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SUMMARY 


A finite element method for the spatial discretization of the dynamic equa- 
tions of equilibrium governing rotary-wing aeroelastlc problems Is presented. 

The equations of motion are nonself-adjoint, nonlinear, and In partial differen- 
tial form. For this class of problems, variational principles are not avail- 
able. >.s, formulation of the finite element equations Is based on weighted 
Galerl'" residuals. This Galerkln finite element method reduces algebraic manip- 
ulative labor significantly, when compared to the application of the global 
Galerkln method to similar problems. However, more computer time Is spent on the 
numerical calculations. 

To illustrate the application of the Galerkin finite element method, the 
coupled flap-lag aeroelastlc stability boundaries of hingeless helicopter rotor 
blades in hover are calculated. The finite element method is used to remove the 
spatial dependence from the equations. The ensuing set of nonlinear, ordinary 
differential equations is linearized about an appropriate nonlinear static equi- 
librium position. The number of nodal degrees of freedom In the discretized sys- 
tem is reduced significantly through a normal mode transformation. The nonlinear 
static equations, determining the equilibrium position, are solved Iteratively 
using the Newton -Raphson method. The linearized dynamic equations are reduced to 
the standard eigenvalue problem from which the aeroelastic stability boundaries 
are obtained. 

The convergence properties of the Galerkin finite element method are studied 
numerically by refining the discretization process. Results Indicate that four 
or five elements suffice to capture the dynamics of the blade with the same accu- 
racy as the global Galerkin method. However, for a reliable analysis, two modes 
for each elastic degree of freedom are required, since the second lag mode deter- 
mines system stability for certain values of elastic coupling. 

Next, the method Is applied to the more practical couplea flap-lag-torsion 
aeroelastic stability and response problem of hingeless helicopter rotor blades 
in trimmed forward flight. Emphasis Is placed on consistent discretization of 
the torsional degree of freedom. 

No previous finite element solutions for the stability and response of non- 
linear, nonconservative systems with periodic coefficients are available. There- 
fore, the general formulation Is specialized to the coupled flap-lag problem In 
forward flight which is used to establish the computational feasiblity of the 
Galerkin finite element method in the forward flight regime. 

The nonlinear, periodic coefficient, finite element equations are linearized 
about a nonlinear time dependent equilibrium position, namely, the steady-state 
response of the system. This response is obtained iteratively using quasilinear- 
ization. Aeroelastic stability is determined from the linearized perturbation 
equations using Floquet theory. 
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SECTION 1 


IHTROIXJCTICN AND REVIEW OF 
PERTINENT LITERATURE 


1.1 Introduction 

Rotary-wing aircraft are widely used today in loth civilian and 
military versions . However, performance characteristics must be 
improved to meet requirements on improved speed, range, payload, 
maneuverability, maintenance <*ad comfort. Of the many components 
governing the performance of rotary- wing aircraft, the rotor is of 
outstanding importance. Consequently, one of the most active research 
areas in aeroelasticity today is rotary- wing aeroeiasticity. 

The coupled flap (bending out-of-plane of rotation), lag (bending 
in-plane of rotation), and torsional aeroelastic behavior of an iso- 
lated rotor blade is the basic building block from which a more com- 
plete system analysis can be developed. A clear understanding of the 
single blade behavior, governed by the complex interaction of struc- 
tural, inertia, and aerodynamic forces, is therefore imperative. Both, 
stability and response, must be well understood. Terminology and con- 
figuration parameters associated with hingeless rotor blades which 
have become an increasingly attractive concept are given in Figs. 1,2. 

Several studies have derived equations capable of simulating the 
motion of this configuration with varying degrees of sophistication. 

A comprehens ive review of recent developments in this area is given by 
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Friedmann [l] . The moat significant conclusion from previous research 
is the fact that the rotary-wing aeroelastic stability problem is 
inherently nonlinear. As a consequence, the correct treatment of this 
aeroelastic stability problem requires the derivation of the dynamic 
equations of equilibrium in a careful and consistent manner such that 
moderate deflections, based upon the assumption of small strains and 
finite slopes, are properly incorporated in the mathematical model (l ] , 
[2 ] . When the equations of motion are formulated in this manner, non- 
linear terms can appear in the structural, inertia and aerodynamic 
operators associated with this aeroelastic problem and the final 
equations of motion will hs'e a partial differential nonlinear form 
[ 1 ], [ 2 ]. 

In rotary-wing aeroelaaticity the nonlinear equations of motion 
in partial differential form are usually solved by applying Galerkin's 
method to eliminate the spatial dependence of the problem [l] - [3]. 

This procedure yields a set of coupled nonlinear ordinary differential 
equations for the dynamics of the blade. It is common practice, [l] - 
[3l, to obtain actual aeroelastic stability boundaries by linearizing 
the equations of motion about an appropriate equilibrium position and 
extracting stability information from the eigendata associated with 
the linearized system. 

Typical studies, [l] - [3], dealing with practical blade configura- 
tions in hover, or in forward flight, are representative of the alge- 
braic complexity encountered when applying Galerkin'3 method to 
rotary-wing aeroelastic problems. From the inspection of these and 
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similar studies it is clear that methods of solution based upon the 
modal Galerkin method lead to extremely cumbersome algebraic manipula- 
tions , which have to be carried out manually or by alternative means 
such as algebraic manipulative systems. In some cases, the amount of 
algebraic manipulations associated with the global Galerkin method is 
go excessive as to prohibit treatment of complicated blade configura- 
tions in a re&listie manner. Therefore, in this study the spatial 
dependence will be eliminated using a Galerkin type finite element 
method. This essentially local Galerkin method enables one to dis- 
cretize the partial differential equations of motion directly. Conse- 
quently, a significant reduction in the algebraic manipulative labor 
required for the solution of the problem is accomplished. 

During the past fifteen years, the finite element method has 
undergone explosive growth and, at the present} it has evolved from a 
structural analysis tool to a general mathematical method for solving 
partial differential equations, which is competitive with finite dif- 
ferences, for general applications, and superior to finite differences 
in structural dynamics applications [4] - [8]. For conservative self- 
adjoint, linear problems, the finite element model for the system can 
be conveniently generated by applying appropriate variational princi- 
ples. Existence of these variational principles will also, in most 
cases, guarantee the com-v.-gence of the method. For nonself-adjoint, 
nonconservative problems, such as the flutter or aeroelastie problem, 
variational principles are not available. Thus, generation of the 
finite element model for aeroelastie, nonconservative systems is more 
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complicated and convergence of the method is not guaranteed 19] • 

The rotary-wing aeroelastic problem is nonself-adj oint , noncon- 
servative and nonlinear, thus, formulation of a finite element method 
for this problem i3 by no means straightforward. However, finite 
eleraen.. discretization of these equations of motion will essentially 
eliminate the cumbersome algebraic manipulations associated with the 
global Galerkln method. 

The purpose of the present study is to develop a local Galerkin 
method of weighted residuals ( 5 ] - ( 8 ], [10], 111], which i3 used to 
discretize the spatial dependence of the equations resulting in a 
finite element formulation of the rotary- wing aeroelastic problem. 

This method is applied directly to the equations of motion in partial 
differential form and leads to a finite element formulation of the 
rotary-wing aeroelastic problem, avoiding the excessive algebraic man- 
ipulations required by the application of Galerkin' s method when using 
global modes (i.e., conventional method). 

To illustrate the method and establish its convergence properties, 
the method is applied to seme typical rotating blade free vibration 
problems and to the coupled flap- lag aeroelastic stability calculation 
of a hingeless helicopter rotor blade in hover. Comparison of the 
solutions obtained, by using the finite element method, with previously 
published results is used to establish the convergence properties of 
the method. It is concluded that this formulation has the potential of 
becoming a powerful and practical tool for solving rotary-wing aero- 
elastic stability or response problems. 
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Next, the Galerkin finite element method, ia extended to dis- 
cretize the coupled flap- lag- tors ion aeroelastic stability and 
response problem of hingeless helicopter rotor blades in forward 
flight. In this case, the equations of motion contain periodic time- 
varying coefficients. For simple blade configurations, finite element 
solutions for the flap- lag case are compared with previously published 
results. In addition, more complicated blade configurations can be 
treated efficiently. Here, the f ull potential of the new method 
becomes evident, as it enables one to model the blade more realistic- 
ally to a degree which has not been achieved previously. 

1.2 Review of Discretization Procedures Used in Rotary-Wing 
Dynamics and Aeroelastic ity 

The rotary-wing aeroelastic problem is governed by partial dif- 
ferential equations. The first step in solving these equations is to 
discretize the spatial dependence of the dependent variables such 
that a set of ordinary differential equations is obtained. Analysis 
of these equations will yield the dynamic system behavior. 

Typically, Galerkin' s method of weighted residuals is used for 
the discretization procedure [l] - [?]. Thus, each elastic degree of 
freedom is represented as a finite sum of mode shapes. These modes 
are taken as the coupled [3] or uncoupled [2] free vibration modes of 
a rotating blade. They arc generated from the uncoupled free vibra- 
tion mode shapes of a nonrotating blade, for which exact expressions 
are available. Most studies, in particular those considering the 
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forward flight case, use only one mode for each elastic degree of 
freedom. When more than one mode is vised, the complexity of the 
problem leads to extremely cumbersome algebraic manipulations which, 
in some cases, sure so excessive as to prohibit realistic treatment of 
complicated blade configurations. 

Thus, it becomes compulsory to look for alternative discretiza- 
tion procedures. Two such methods, the integrating matrix method and 
the finite element method, will be introduced in this section. There 
are, of course, other methods. However, these two have been applied 
more widely to rotating beam problems and seem to have greater promise 
with respect to treatment of aeroelastic problems. 

The integrating matrix method (IMM) provides a means to eliminate 
the spatial dependence in rotating beam vibration problems. Vakhitov 
developed the method and used it to solve for static deflections of 
beams [12] and for coupled bending- tors ion vibrations of a rotating 
blade [123 • However, only few numerical results were presented. 

Hunter [ lh ] extended the numerical scope of the method and investigat- 
ed coupled bending-bending vibrations of a rotating propeller blade. 
Subsequently, several researchers applied the IMM to a variety of 
rotating beam vibration problems. White [15] formulated the coupled 
bending- tor' ion problem. Murtby [l6] investigated flapwise bending 
and coupled flap-torsion vibrations. White and Malatino [l8] solved 
the flap-lag and the nonlinear torsion problem for the same propeller 
blade as considered by Hunter [lh] . Finally, Kvatcrnik, White, and 
Kaza [19] used the IMM to solve nonlinear flap-lag and axial vibration 
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problems. They also indicated the possibility of applying the IMM in 
aeroelastic stability analyses of helicopter rotor blades. In what 
follows, a brief description of the H44, its basic properties and 
numerical performance when applied to rotating beam problems, is 
given. 

The IMM is based on direct numerical integration. The integrating 
matrix may be viewed as a matrix operator which by premnltiplying a 
vector, containing as elements the values of a function at discrete 
stations along the blade, transforms it into another vector having the 
integrals of the function (from one end of the blade to each station) 
as elements. To account for the boundary condition, a constant 
vector has to be added. In order to apply the IMM towards the solu- 
tion of a differential equation, it is necessary to write the differ- 
ential equation, or an integrated form of it, at a number of stations 
along the blade. The resulting set of equations has to be cast in 
matrix form. The integrating matrix can then be used to express the 
equations in terms of one set of unknowns, either the displacements or 
the fundamental derivatives (second-order for bending, first-order for 
torsion) at each station. Thus, discretization is achieved and the 
vibration problem is now posed in the form of a matrix eigenvalue 
problem. 

Derivation of the integrating matrix is based on piecewise poly- 
nomial interpolation. If, for convenience, equally spaced collocation 
points are chosen, Newton's forward-difference interpolation formula 
can be used to express the polynomial coefficients in terms of the 
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function values at the appropriate collocation points. Integration 
of the polynomial expressions yields the elements of the integrating 
matrix. 

The IMM always leads to nonsynsaetric system matrices, even when 
considering a self-adjoint problem. Furthermore, the matrices are not 
banded. The eigenvectors are not orthogonal with respect to the sys- 
tem matrices. Further, the IMM does not yield upper bound solutions. 
Last, the dynamic matrix is degenerate, leading to zero eigenvalues 
which correspond to infinite frequencies. 

The inputs for the matrix equations are simply the values of the 
cross-sectional properties at the discrete stations. Nonuniform 
properties axe therefore easily incorporated. Boundary conditions are 
readily applied when considering the clamped- free case. For other 
cases, modification of the method is necessary [12 3 and some of its 
appealing simplicity i3 lost. 

The IMM has been applied to a number of static, vibration, and 
buckling problems of beams. Results were compared with exact and 
other approximate solutions, employing the finite difference, transfer 
matrix, and Rayleigh- Ritz methods, and with experimental results. 
Extensive convergence studies were performed by Hunter [ih-] for bend- 
ing free vibrations of a cantilever beam. Overall, accuracy increases 
with the number of stations and the degree of the interpolating poly- 
nomials. In general, higher degree polynomial representation gives 
higher accuracy when using a fixed number of stations. However, the 
number of stations employed should always oe considerably larger 
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(afccut twice as much) a. the degree of tha interpolate polynomials, 
especially when higher frequencies are desired. A convergence study 

f°r " , ° re 0Onplex 3 >' 3tK “' “Wch for instance display dynamic in- 
stability, was not performed. Results for nonlinear problems are 
indicative of trends only. 

lUe DM has been established as a very useful tool to solve free 
vibration problems of rotating beams. The method is very accurate and 
converges well. This is not surprising as it is solely based on 
numerical integration. The Mi has the potential of being applied to 
rotary-wing aeroelastlc problems. However, more information on its 
performance, when solving nonlinear and nonconservative problems, 
needs to be obtained. One particular disadvantage of the method is, 
that it does not yield orthogonal mode shapes. Thus, it becomes 
questionable whether modes free, the MM could be used to reduce the 

number of degrees of freedem, which is essential for calculating aero- 
elastic stability boundaries. 

”” — eleEent (f™) has been used extensively and 

With great success for the solution of been vibration problems. In 

particular, it has been applied in the analysis of rotating blades. 
Solutions for the flapwiee free vibrations of a pinned-free rotating 
beam using a Galerbln-type finite element method were presented by 
Nagaraj and Shantahumar [20], Their treatment was based on a bending 
element with four degrees of freedom per node, satisfying all boundary 
conditions. Comparison with results obtained by using the global 
GalerMin method showed good agreement. A conventional finite element 
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model for the natural vibrations of tapered and pretwisted canti- 
levered rotor blades was formulated in Reference [2l] . The equations 
were developed for coupled flap- lag- tors ion motion. A simple finite 
element, having ten degrees of freedom, was employed. However, numer- 
ical results were restricted to flapwise bending of a rotating tapered 
beam and a nonrotating pretwisted beam. The results displayed satis- 
factory convergence trends and agreed well with previous solutions. 

Application of the FEM to nonlinear problems is, by now, 
standard. Many solution algorithms are available, and their numerical 
performance is well documented. In addition, the FEM has been used 
for the analysis of a large variety of nonconservative systems. Here, 
the Galerkin finite element method has found its most prominent area 
of application. 

The FEM, in contrast to the IMM, can always be formulated such 
that it yields symmetric matrices for self-adjoint problems. This, 
together with the banded nature of the matrices represents an advan- 
tage for certain solution algorithms. Another advantage of the FEM is 
the handling of boundary conditions of any type without any modifica- 
tions. On the other hand, formulation of nonlinear terms and, in 
particular, integral terms, as they appear in rotor blade equations, 
is quite straightforward when using the HIM. Finally, the FEM leads 
to orthogonal eigenvectors. Such eigenvectors have been successfully 
employed in coordinate transformations to reduce the number of degrees 
of freedom used to model a dynamic system. Such a reduction is of 
great importance, since the aeroelastic problem at hand requires a 
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large number of eigenproblem solutions. This last point and the 
large amount of information available on application of the FEM to 
nonlinear and nonconserv&tive systems, represent the main advantages 
of the FEM over the HIM. It also should be noted that the IMM is 
purely a nume r ical integration procedure, while the FEM retains a re- 
lation to the physical properties of the system under consideration. 

A thorough examination of the relevant FEM literature will fol- 
low in the next section. From this, it becomes clear that the Galer- 
kin finite element method is extremely well suited for dealing with 
rotary-wing aeroelaaticity. A further advantage in using a Galerkin- 
type finite element approach consists of the considerable amount of 
research done by applied mathematicians and engineers to establish the 
numerical properties of this method. This vigorous, ongoing research 
activity provides the aeroelast ic ian with more information on the num- 
erical aspects, and particularly, convergence properties of the method 
than is available on the integrating matrix method. 

1.3 Review of I^rtlnent Literature on Finite Elements 

From previous remarks, it is clear that two aspects of the finite 
element method are of particular concern here. Namely, application to 
problems which do not allow a variational formulation and nonlinear 
problems. Accordingly, the following literature review will emphasize 
these special features. 

The finite element method (FEM) originated in the area of 
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structural analysis. Its formulation was based on elementary concepts 
such as the principle of virtual work or Castigliano's theorem. 
Recognizing that the FEM can be viewed as an application of the princi- 
ple of stationary potential energy soon led to its use in other areas. 
In general, all problems for which a variational principle existed 
could be solved. This variational formulation also, in moat cases, 
guarantees the convergence of the method with mesh refinement. How- 
ever, there is a large number of problems for which no variational 
principle exists. 

Olson [22] solved the nonself -adjoint panel flutter problem. He 
used the principle of virtual work to derive consistent aerodynamic 
load matrices. The eigenvalues did converge when the number of ele- 
ments was increased. However, this convergence was not monotonic. 
Barsoum [2J] used the extended Hamilton's principle to solve the dyna^ 
mic stability of thin-walled beams subjected to nanconservative 
loading. The nonself -adj oint character of the system was reflected in 
the nonsymmetric load matrix. The same problem was treated by 
Kikuchi (9]» It is interesting to note that convergence of tne re- 
sults was best when the system was conservative, but became worse as 
the degree of nonconservativeness became dominant. 

A much more general approach was first developed by Szabo and Lee 
[2L]. They combined the method of weighted residuals, using Galer- 
kin’s weighting criteria, with the FEM to calculate stiffness matrices 
for problems in plane elasticity. This procedure requires only know- 
ledge of the differential equations and boundary conditions in a given 


12 



domain and boundary. No variational concepts are involved. In Refer- 
ence [2h] the weighted residual of the governing differential equa- 
tions was evaluated over the element (local domain) . Integration by- 
parts led to a formulation which yields the same element matrices as 
the conventional FEM. This also introduced element boundary terms, 
which .would vanish in the element assembly process if (and only if) 
both the approximate displacement and stress fields were continuous. 
This is, in general- not the case. Therefore, inter-element boundary 
contributions were neglected while, for external boundaries, the act- 
ual boundary conditions were substituted. 

Zienkiewicz and Rarekh [25] applied the Gale r kin weighted resid- 
ual finite element method (GFEM) to transient field problems. In 
contrast to Reference [2k], the residual was formulated for the system 
(global domain). Green's theorem was used to remove higher order con- 
tinuity conditions between elements. Formally, no inter-element 
boundary terms appeared because the integration was carried out over 
the global domain. 

A more rigorous treatment of the GFEM was presented by Hutton and 
Anderson [ll] . They used approximating functions over the global do- 
main which were nonzero only within the local domain. This made it 
possible to apply convergence results from the global Galerkin method. 
Further, it was clearly stated which of the boundary conditions have 
to be satisfied and what the inter-element continuity requirements 
are. Both issues are closely related to the integration by parts. It 
was also shown that inter-element boundary residuals need to be 
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Introduced when deriving the equations on the element level. Finally, 
it was established that all those problems that can be solved using a 
variational approach are a subclass of those amenable to the GFEM. 

If a variational principle is available, both methods lead to the 
same results. Last, deriving the equations on the basis of virtual 
work is equivalent to using Galerkin's method as long as the geometric 
boundary conditions are homogeneous. 

Aral, Mayer and Smith [26] discussed the relation between G&ler- 
kin's approach and true- , quasi- , and restricted variational princi- 
ples. If available, they all lead to the 3ame results. The approxi- 
mating equations, using the GFEM, were derived on the element level. 
Consequently, inter-element boundary residuals were introduced to can- 
cel identical terms arising from integration by part 3 ever each 
element . 

Use of the GFEM without performing integration by parts was il- 
lustrated in References [27] - [29]. Prasad and Murty [27] solved 
flexural beam vibration problems. They used seventh-order interpolat- 
ing polynomials as shape functions to satisfy continuity requirements 
on the higher derivatives. Comparison with the conventional FEM (cub- 
ic interpolating polynomials) proved the GFEM, in this fora, to be 
very accurate. In Reference [28] the stability of nonconservative 
systems was analyzed. Results were more accurate than those of 
Kikuchi [9]* This improved accuracy has to be attributed to the stip- 
ulation of higher-order continuity, see Barsoum [23]* A disadvantage 
of this particular implementation of the GFEM becomes apparent when 



considering problems with mixed boundary conditions. Then it is 
necessary to perform a congruent coordinate transformation such that 
these boundary conditions can be satisfied [28] , [29l. Furthermore, 
matrices which would be symmetric when applying integration by parts 
now become unsymmetric. 

A comparative study of several finite element models applied to 
vibrations of beams was presented in Reference [30] . The basic fea- 
tures of the conventional- , hybrid- , least-square- , collocation- , 
and GFEM were discussed. The possibility of using the GFEM in non- 
linear analyses was pointed out. Further material on several approx- 
imation procedures, used in conjunction with the FE24, can be found in 
Reference [8]. 

Application of the FEM to nonl-inear structural problems has re- 
ceived considerable attention and a large number of publications on 
this subject are available. A comprehensive review of solution pro- 
cedures applied in static analyses of structures, displaying both 
geometric and/or material nonlinearities, was presented by Tillerson, 
Stricklin and Haisler [31] - Selection of a solution procedure was 
shown to be governed by interaction of several factors, such as type 
of analysis, ease of implementation, storage space, problem size, 
desired accuracy, degree of nonlinearity, computational economy, and 
user experience. Another review by Gallagher [32] dealt with geometri- 
cally nonlinear problems. Construction of finite element equations 
and solution algorithms were seated. It was shown that tensor nota- 
tion [33], rather than matrix notation, permits a greater simplicity 
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and efficiency in formulating the equations. It was also pointed out 
that some advantage might be gained by using inconsistent formula- 
tions, in which simplified interpolation polynomials are used for com- 
putation of the nonlinear terms. Such an approach was used by Bergan 
and Clough [34]. A more detailed investigation of inconsistent formu- 
lations can be found in Reference [35] • Finally, Gallagher [32] sug- 
gested use of condensation techniques to reduce the number of degrees 
of freedom, prior to performing the nonlinear analysis computations. 
This could lead to significant savings in computational expense. The 
normal mode concept was applied for this purpose by Kavanagh [36]. 

The FEM has also been used successfully in problems of nonlinear 
structural dynamics. However, applications were mainly restricted to 
the transient response tinder impulsive loadings, such as impact, 
seismic, or blast loads, wherein the time history was found by numeri- 
cal integration [37] - [393 * 

The first attempt to solve large amplitude natural vibrations of 
beams and plates was made by Mei [4o], [4l]. He calculated the non- 
linear terms from the linear mode shape, multiplied by an amplitude 
factor, and then extracted the eigendata from the linearized system. 
Comparison with other analytical solutions and experimental data 
showed the FEM to match the experiments more closely. Convergence is 
mcco ionic with increasing mesh refinement. Mei, later on, included an 
iterative solution technique, where the nonlinear terns were calcu- 
lated from the vibration mode shape of the previous iteration cycle 
[ 42 ], [43]. The iteration led to lower frequencies. The number of 
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iteration steps to achieve a certain accuracy was larger for higher 
values of the prescribed amplitude, i.e., it increased with the degree 
of nonlinearity. Convergence was not monotor.ic during the iteration. 
The same problems were solved in References and [I 5 ] using the 

same iterative procedure. However, the nonlinear terms were modeled 
in a different manner, without making some of the simplifying assump- 
tions used by Mei. Comparison of the period for the fundamental mode 
of a 3imply supported square plate showed results of Reference [i-5] to 
be more accurate than those of Reference [ll]. This illustrates the 
importance of modeling the nonlinear terms adequately. 

Of particular interest is Mel's finite element approach to non- 
linear panel flutter r 16 ] which represents a typical nonself-ocijoint 
nonlinear aeroelastic problem. Development of the aerodynamic matrices 
was based on Olson's work [22]. Modeling of the nonlinear terms, the 
iterative solution procedure and equivalent linearization technique 
were taken from References [^l] and ( U2 ] . For large deflections, the 
nonlinear effects, mainly due to membrane stresses, restrain the panel 
motion to bounded limit cycle oscillations with increasing amplitude 
as the dynamic pressure increases. In general, three to six iteration 
cycles were sufficient for convergence. Convergence with respect to 
the number of elements used was not studied for the nonlinear case. 

The FEM described the panel behavior correctly. However, a numerical 
comparison with previous nonlinear panel flutter solutions was not 
given. Thus, the accuracy of the nonlinear finite element solution 
cannot be assessed. 
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Finally, it is worthwhile to mention that the GFEM is very well 
suited for application to nonlinear problems. Oden [17] characterized 
it as "... perhaps the most powerful technique for generating accept- 
able finite element models of nonlinear equations." Rao and Raju [kfl] 
applied it to the post buckling analysis of uniform cantilever columns. 
Agreement of results with those found by using elliptic integrals was 
very good, even for a small number of elements. The GFEM was also 
used in nonlinear reactor dynamics [k9l and nonlinear boundary layer 
flows [50]. This illustrates, again, the versatility of the GFEM. 

l.k Objectives of the Present Study 

The objectives of the present study are summarized below: 

1. Development of a local method of weighted residuals, using 
Galerkin'a weighting criteria. This results in a finite 
element formulation which can be used to discretize nonself- 
adjoint and nonlinear partial differential equations 
encountered in rotary-wing aeroelasticity. 

2 . Application of this Galerkin finite element method (GFEM) to 
coupled flap- lag aeroelastic stability calculati . of hinge- 
less helicopter rotor blades in hover. The GFEM is used to 
transform the nonself-adjoint, nonlinear partial differential 
equations of motion into a set of ordinary nonlinear differ- 
ential equations. The equations are then linearized about 

an appropriate nonlinear static equilibrium position. The 
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number of unknowns modeling the discretized system is re- 
duced significantly through modal analysis. The nonlinear 
static equations are solved iteratively, using the Newton- 
Raphson method. The linearized dynamic equations yield a 
p isuulard eigenvalue problem from which the aeroelastic 
stability boundaries are obtained. The convergence proper- 
ties of the GFEK are studied numerically by refining the 
disc: ' ,J .zation procedure. Results will also be compared 
with previous analyses where the global Galerkin method was 
employed. It should be stressed here that the major inter- 
est i 3 the spacewise discretization of a typical rotary- 
wing aeroelastic problem via the GFEM. 

3. The Galerkin finite element method is used to discretize 
the coupled flap-lag- tors ion aeroelastic stability and re- 
sponse problem of hingeless helicopter blades in trimmed 
forward flight. Emphasis is placed on consistent discreti- 
zation of the torsional degree of freedom. The blade is 
assumed to have built-in twist, cross-sectional offsets be- 
tween the aeroelastic center, center of gravity, and elastic 
center, and nonuniform mass and stiffness properties. Root- 
torsional spring stiffness and nonuniform cyclic inflow are 
also included in the model. When solving the discretized 
dynamic equations, a major complication arises due to the 
forward flight condition which introduces periodic time- 
varying coefficients in the equations of motion, as well as 
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a large number of additional aerodynamic loading terms. 

Due to the complexity of this problem it was decided, 
at this stage, to specialize the general formulation to the 
coupled flap-lag problem in trimmed forward flight. This 
problem, which is computationally less expensive, is used to 
establish the computational feasibility of the Galerkin finite 
element method for the aeroelastic problem in forward flight. 
m this portion of the present study, the discretized, non- 
linear ordinary differential equations of motion are linear- 
ized about a nonlinear time-dependent equilibrium position, 
namely, the steady-state response of the system. Aeroelastic 
stability boundaries are obtained from the linearized system. 

For a number of sample problems, involving hingeless 
rotor blades, the finite element solutions are compared with 
previously published results obtained by usir,g the global 
ualerkin method. The convergence properties of the Galerkin 
finite element method are studied numerically by varying the 
number of elements and mode shapes, i.e., by refining the 
discretization procedure. The sensitivity of the aeroelastic 
steady-state response and stability to variations in the 
parameters governing this problem is also considered. 
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SECTION 2 


GALERKH; FINITE ELEMENT METHOD 


Froo the preceding literature review it is evident that the 
Galerkin finite element method (GFEM) has been used successfully in 
treating a large number of nonself -adjoint and nonlinear problems. In 
the next section a brief description of the method will be given. 
Emphasis is placed on its application to nonself -adjoint systems. The 
GFEM, being applied directly to the governing differential equations, 
makes discretization of nonlinear terms straightforward. The ensuing 
nonlinear equations can then be solved with any of the algorithms used 
in the conventional finite element method. This aspect is treated in 
Sections 3*3 and b. 3 . 1 . 

2,1 Global Galerkin Method 


The local Galerkin method, resulting in a finite element discre- 
tization, can best be clarified by illustrating its application to a 
simple system, f-lore details can be found in References [8] and [UJ . 

Consider the following differential equation 

Q(q) + P(q) = F (2.l) 

which is defined in a domain D, where Q is a symmetric positive 
definite differential operator of order 2r and P is a general 
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is assumed. The • are linearly independent shape functions and the 

b are the undetermined parameters for this problem. When using the 
m 

extended Galerkin method, the * have to have continuous deriva- 

m 

tives up to order (r-l), i.e., C r _^ continuity. Further, they 

need to satisfy only the geometric boundary conditions, i.e., those 
containing derivatives of order not higher than (r-l). This approx- 
imate solution is then substituted into the differential equation and 
the boundary conditions. The error is minimized by requiring ortho- 
gonality with respect to a set of weighting functions. Thus, an 
integral statement, equivalent to the differential equation and the 
boundary conditions, is obtained. 

In the extended Galerkin method the original shape functions 
are chosen as weighting functions. It is then required tint the sum 
of the weighted residuals of both the differential equation and the 
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natural boundary conditions, i.e., those containing derivatives of 
order r and higher, be zero [10] . Thus, 


♦ 6 dD + 
m 



/ * e_dS»0, m = 1,2 , . . . ,M , (2.3) 

J o mb 


where 

£ = Q(q g ) + ?(q g ) - F (2.4) 

and £_ is the residual associated with the natural boundary 

B 

conditions. 

In many cases it is possible to apply integration by parts to 
Eq, (2.3). This reduces the order of differentiation in the symmetric 
operator Q, thus lowering requirements on the shape functions* 

m 

from CL to C continuity. Furthermore, this algebraic step 
also yields terms which cancel some of the boundary residual contri- 
butions. As a matter of fact, formulation of the weighted boundary 
residual and its combination with the differential equation residual 
are made in such a way that, when integrating the last one by parts, 
identical terms of the boundary residual are canceled [10]. In the 
case where the natural boundary conditions are homogeneous, all 
boundary residual terms cancel out and Eq. (2.3) becomes 

J {Q(* m ,q g ) + * m F(q 8 ) - \ F}dD = 0 , (2.5) 

where Q denotes the operator Q after integration by parts. 
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Application of the extended Galerhln method to beam bending can be 
found In Reference [ 51 ] . 


^ Finite Element Approach 

When formulating a finite element version of Galerkin's method 
the danain D is subdivided into E subdomains d, which are 
called elements. In each element an approximate solution of the form 
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is assumed, where a* are the nodal parameters and *• are linearly 
independent shape functions defined only in the subdomain associated 

with the element d. This local approximation can be extended over 
the whole domain D by defining 
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Using Eq. (2.7), the global approximation can be expressed as 
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After imposing compatibility conditions on the nodal parameters 
of adjacent elements (this is done during the process of assembly), 
equations (2.2) and (2.8) are equivalent. 

Equation (2.5) can be rewritten using Eq. (2.8) as 

X X <5 - S F j a °* <2 - 9) 

i=>l n*l ' 

J 3 1,2 , . . . ,N, e 3 1, • • • ,E, 

where 3. is obtained from Q by means of the previously mentioned 
integration by parts. In addition, it is implicitly assumed that no 
inter-element discontinuities occur. Thus, the functions and 

its derivatives of order less than r must be continuous such that 
derivatives from order r up to (2 r-l) are finite on element 
interfaces. Equation (2-9) represents all weighted residuals for one 
element and for the total assemblage of elements (i.e., N*E 
weighted residuals), it represents an intermediate step and in reality 
only Eq. (2.10) below is used. 

As a consequence of the linear independence of the ^ functions. 
Equation (2.9) can also be rewritten on the element level. 



J 3 1,2 , . . . ,N, e 3 1,2 , . . . ,E . 
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Equation (2.10) thus represents a set of N equations for each 
element, from which the element matrices can be calculated. 

Because of the unique choice of weighting functions in Galerkin's 
method and because of the integration by parts, Q yields symmetric 
matrices. Further, when P is equal to zero, Eq. (2.10) is the 
exact same expression as found when employing the variational formu- 
lation of the finite element method. Operator P leads to unsym- 
metric matrices. However, the banded nature of the system matrices is 
still preserved. 

Assembly of the system matrices and enforcement of the geometric 
boundary conditions is handled as in the conventional finite element 
method. 

F in ally, it should be pointed out that when Equation (2.5) is 
solved directly, the approximate solution has to have C, r contin- 
uity and must satisfy all boundary conditions. The generation of 
such finite elements is of course more difficult, in particular for 
nonlinear terms, than generation of elements for the solution of 
Equation (2.5). In addition, all matrices will be nonsymmetric . On 
the other hand, due to the higher-order continuity, one might expect 
more rapid convergence. Thus, it becomes obvious that integration by 
parts plays an important role. 

2 . 5 Convergence Properties 

The Galerkin finite element method is equivalent to the 
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conventional finite element method when considering self-adjoint 
problems. It is well known that elements which are conforming and 
are able to approximate constant strain will ensure convergence for 
this class of problems. Some elements even display monotonic converg- 
ence, thus allowing use of efficient extrapolation procedures and 
give an upper bound on the potential energy. 

Based on Mikhlin's work [ 52 ], [ 53 ], Hutton and Anderson [ll] and 
Kikuchi [ 9 ] established convergence criteria for the Galerkin finite 
element method when applied to a wider class of problems than those 
amenable to the variational FEM. However, numerical results show that 
convergence is, in general, not monotonic [22] and becomes less rapid 
when the nonself-adjoint character of the system under consideration 
becomes more pronounced [ 9 ]. 

Convergence studies for the Galerkin finite element method, when 
applied to nonlinear systems, are of numerical nature only [48j. 

Noor and Whiteman [54] derived an error bound for a certain class of 
nonlinear problems, solvable by the GFEM. There are a number of 
studies on convergence, accuracy, and stability of the FEU in non- 
linear problems. They are, however, either too general to be useful 
for practical applications or restricted to certain special classes of 
problems. On the other hand, a large number of nonlinear problems 
have been solved using the FEM with great success. 

The comments made in this section are basically to be understood 
as an indication of the ongoing research effort. The rotary-wing 
aeroelastic problem — because of its complexity — will hardly be 
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accessable to any convergence proof. Thus, for the time being, 
convergence can only be established numerically [55] > i.e., by refin- 
ing the discretization process. 
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SECTION 3 


APPLICATION OF THE METHOD TO THE 
FLAP-LAG AERGELASTIC PROBLEM 
IN HOVER 


3.1 Brief Description of the Coupled Flap- Lag Equations 

of Motion 

The coupled flap-lag equations of motion used in this study- 
serve mainly as an Illustrative example for the application of a 
Galerkin-type finite element method to rotary-wing aeroelasticity. 

The flap- lag equations for hover are obtained from the general equa- 
tions which have been presented in Reference [2 ] , by an appropriate 
elimination of the terms associated with forward flight and the 
torsional degree of freedom. 

The geometry of the problem is shown in Figures 3 and 4. A few 
important assumptions made in the derivation of these equations are 
briefly stated below: 

1. The blade is assumed to have moderate deflections, which 

implies small strains and finite rotations or slopes. These 
elastic rotations are assumed to be of order (e^ - 0.20) 

30 that terms of 0( e D ) are negligible compared to terms of 
order one , 0 (l ) . The blade can bend in two mutually perpen- 

dicular directions. Initially the blade is straight; 
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during deformation the Euler -Bernoulli assumption is used. 
The structural operators resulting from these assumptions 
have been presented in Reference [56]. 

2. The blade has only precone (5^; it is cantilevered to the 
hub and there is no built-in twist. Aerodynamic, tension, 
mass, and elastic centers coincide, i.e., all cross- 
sectional offsets are zero. Inertia and stiffness proper- 
ties are assumed to be uniform. 

3* There is no coupling between blade and fuselage dy nami cs, 
i.e., pitch, roll, and yaw motions shown in Figure 1 are not 
coupled with the dynamics of the blade. 

4. Two-dimensional quasi3teady aerodynamic loads are used; 
apparent mass, stall, and compressibility are neglected. 

5* An ordering scheme identical to the one in Reference [2] is 
used and quantities having the magnitude of the squares of 
the blade slopes are neglected when compared to one, i.e., 

0(1) + 0(e Z ) - o(l) . 

This ordering scheme is given in Section 4.1. Note, that in 
Section l>, e^ = O(e^) is assumed. 

Using these assumptions, the coupled flap- lag equations of motion 
for hover can be written as follows . 

Axial equilibrium : 

* 

*, x + <V« X ) (}.l) 
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Flap equilibrium : 


- li,, V - B,, w + (w - 6 CH + e ) - 28 v 

6 3 j XX30C 33 > xxxx ,x ,x P p vx 0 e l y ^ p p v 

♦* 

- W + + 2 « 1 )- (x 0 + «i) 7 ^ - Xq p p v + (^ ex o"7 ^)v 

* 

" ( *0 + *1 } * - *0 * ", x + *0 7/ ,x \x ' ^ " 0 (3.2 ) 


La^ equilibrium : 


^ 2 V ^xxxx ^3 W ,xxct + (v ,x ^,x + d* J (v , x + w ,x )d ^0 


+ v-v + 2 Bp v. r[-s(^. ^J-x* ( £ xr-JS + 2 ^) 

* (9 i x * 

+ [2 7 X - ©(j^ + e 1 )]w + (2 7 X - Gx^v w ^ 


x jL v 


x n v w +2£J vw-Gvw+v 

u ) A |X U 


-2 


+ 2vw W-X-.V w wl 
>x V ,X ,x J 


(3.3) 


at 


at 


The boundary conditions are given by [56] 


*0 * 0; 


V*W=»V a w a 0 , 

>X ,x ’ 


Xq 3 1; 


Bp., V + EL, v a 0 
22 ,XXX 23 , XXX 


(3-W 


(3* 5a) 
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Bn, 

V + 

B„ 

w "0 

(3- 5b) 

^3 

, XXX 

33 

, XXX 


®22 

V + 

,xx 


*:« 

■w 

e 

u 

O 

(3.5c) 


V + 

,xx 

B 33 

o 

n 

B 

•a 

• * 

(3*5d) 

f - 

0. 



(3-5e) 


3.2 Implementation of the Galerkin Finite Element Method 

The Galerkin finite element method, as developed in Section ' 2 , 
is now applied to Equations (3*2) c-nd. (3*3)- To facilitate manipula- 
tion of these equations they are rewritten in matrix operator form. 

- tC^, ^)](4> + ([Gl + t»i^o )] 

+ [^(q, x 0 )]){4} - 

+ ([S B ] - [S T1 (x 0 )] - [K^] + 

+ [^(a, ^)]fq] - fF(x D )^ (3.6) 

where fq] » ( J | . Other quantities are defined in Appendix D . It 
should he noted that th<- tension T was eliminated using the axial 
equat ion . 

According to Eq. ( 2 . 2 ), an approximate global solution is given 

by 
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V ts B>^ - d,]R g } - (K,]£ t *) 

^k. 

» [B] D x {q 8 }-D x ((T 1 ]{q 6 } + [£j{q g }) + - [K^fq 8 } (3.9a) 

Ep " <«*,] + tG] + [D x ] + 0^3 + [D 3 ]){q S } + ([A 1 ] + [Ag] ){*«} (3.^) 

[ V T ~B = < ” C *J T ^ D x^ 8 } + DJM T [B] D 2 {q g } 

X 8k X 

+ fV I([T iHi g ]+ f^)K g }))l (3.9c) 

''o ' 1 

and D* - £/$£ . 

Integrating the term by parts, Eq. (3.8) becomes 
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(3-10) 



l*V 

tn 


(e - f) 

~P 




0 


t 


M x 1 


where 


k 3 #V* [Bl 4t» 8 ) + y*,! 1 dpfq 8 } 

♦yy* !%]«?] + ti.iu 8 ) . t^tA. (3.U) 


Inside the element the displacements sure given by 





(3.12) 


In the present analysis, ordinary beam- type bending [5] . [7] 
representation was chosen for modeling the flap-lag motion; the geo- 
metry of the element is illustrated in Figure 3b , where & e and h e 
represent symbolically the four nodal degrees of freedom associated 
with flapwise and lagwise bending, respectively. Thus, H = k, 

7 n ' \ are Hermite interpolation polynomials. Appendix C, and 

h n and represent nodal displacements and slopes for lag 
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(h n ) » respectively. For convenience, the super- 

script e is omitted from the shape function matrix m. 

As has been indicated previously, stability boundaries in 
rotary-wing aeroelasticity can be obtained by linearizing the equa- 
tions of motion about an appropriate equilibrium position. For the 
case of coupled flap-lag in hover, the equilibrium position is taken 
as the static nonlinear equilibrium position and the dynamic equa- 
tions of motion are perturbed about this equilibrium position, i.e., 

ta e (*)) - {a 50 ] + [i a e (*)] . (3.13) 

Equation (3.12) is now extended over the global domain, in the 
sense of Eq. (2.7) and substituted, together with Eq. (3.13), into 
Eq. (3.10). This yields the nonlinear static equilibrium position 

([IT] + [T*] - [K^l + [aJ] + [A^] ){a e0 } . {F e } (3.1k) 

where it should be noted that fA^] also depends on fa® 0 "} and thus 
Eq. (3.1k) is nonlinear. 

S im i lar ly, the linearized dynamic equilibrium equations are 
given by 
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* Ug'j * f D «) * (I g, „ [d . !)(a ., 5 + 


[x|] {Aa e } 


A 


L [J 2 i] ^ - ^L]{^a e ] 

i-efl Ax 


e-1 


-y 


* <tB e J * f T »] 

i«l 1 


• OP 


+ fAj] + [^])fAa e ] =, {o}. 


6 " 1 ' 2, ***' E • ( 3 . 15 ) 


Equations ( 3 *l 4 ) and (1 « 

, , e _ ^ m equatl °“ Wlttm at the element 

level, thu3 [B e j,fi e ] . 

XJ-.., etc. , represent element matrices which 
ere defined in Appendix D Fo-r +u 

ppe ^d. For the cubic shape functions which have 
B«n Elected, . 

yields the exact element 

Evaluation o f the constant term {p e } ^ ^ ^ 

' *5 * % Md ^ ™rlum equations, 

%1> ' 1S Stra ^ ht “- ^ ™<= term UJI depends on 

tiC eqUUlbrl ““ ~ 1W “ — • Bince the ^ ^ 

position equation will be solved by iteration, [>£] couId be 

evaluated using the value of fa® 0 ] ^ v 

1 j fr ° a tae Previous iterative step. 

ere, . deferent approach is tahen. Rather than calc^atin* (A ?) 
the nonlinear system matrix, Eq. (3.18), la orated directly; see 
F^ndix E . AH the terms in the dynamic equation, Eq. (3.15), are 

[c .~ the “ 3 ^ ^ ^ ^ 

At *2 depend on the static equilibrium position, which is 
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a known quantity after the static problem has been solved. Thus, 
after numerically integrating the triple (quadruple, for [D®]) 
products of the shape functions, again using six-point Gaussian quad- 
rature, the known values of {a 6 ^} are used to evaluate those 
element matrices. 

The assembly of the element matrices into the complete system 
matrices is the same as in the conventional finite element method. 

It should be noted, however, that bandedness of the velocity propor- 
tional matrix is destroyed due to terms of type 

£ and ^ 

i - e+-l i»l 

which are a consequence of the inner integrals over the global domain 

in the operators [C„] and [C, ]. 

Ax 

Here, it should be mentioned that the nonlinear terms were 
modeled completely consistent with the linear terms. No simplifying 
assumptions were involved. This, of course, implies added computer 
time and storage. 

Equations (3-lk) and (3*15) can now formally be written on the 
system level. Thus, the partial differential equations, Eqs. ( 3 . 2 ) 
and (3-3), are reduced to a system of ordinary differential equations 
by using the Galerkin finite element method. 
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3*3 Method of Solution 


Before actually solving the finite element equations, the con- 
siderable number of unknowns representing the nodal degrees of freedom 
will be reduced by applying modal analysis. The basic assumption is 
that the nodal displacement vector can be expressed in terms of a 
small number of mode shapes which approximate the free vibration 
mode shapes of the system. Thus, 


{a 0 } - [A]{q°] 

[Aa] = [A] [Aq] 


(3.16) 


where [A] is the modal transformation matrix containing as columns 
the first M approximate free vibration mode shapes of a rotating 
blade. In the present study, unless otherwise noted, the M^, lowest 
uncoupled nap modes and the lowest uncoupled lag modes are used, 

i.e., M “ + Mp. These mode shapes are also determined by using 

the finite element approach as applied to the rotating beam free 
vibration problem. Since the free vibration equations are self- 
adjoint, the Galerkin finite element method for this case is identical 
to a conventional finite element method. Furthermore, {q 0 } and 
{Aq} are the reduced vectors of generalized coordinates. For the 
static equilibrium position 



JL> 


* J 


0 ,T 


(3-17) 
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and similarly, {Aq] is the reduced perturbation vector of the gen- 
eralized coordinates. The h^,h^,..., etc., can be interpreted as 
amplitudes of the corresponding mode shapes. They constitute the new 
unknowns of the problem. As a general rule, use of M assumed modes 
•will yield m/ 2 actual mode shapes and frequencies with good 
accuracy. 

Clearly, modal analysis provides an effective reduction in the 
size of the eigenvalue problem required for the solution of the 
dynamic system equations. This is a considerable advantage since 
determination of stability boundaries requires repeated solutions of 
the eigenvalue problem. Furthermore, due to this approach, the band- 
edness of the finite element system matrices becomes irrelevant since 
the reduced system matrices are fully populated anyway. Also, the 
geometric boundary conditions can be enforced implicitly through the 
free vibration mode shapes. Finally, it is important to realize that 
modal analysis facilitates the solution of the nonlinear static equi- 
librium equations. Thus, the reduced number of unknowns allows one 
to calculate the derivatives of the nonlinear terms conveniently and 
a relatively efficient solution algorithm, based on the Newton-Raphson 
technique can be used; see Appendix E . 

The final equations, after modal reduction, for static equilir. 
brium are 


[s L ]fq 0 } + C s IJL (q°)Hq°} 11 {cl 


(3.18) 
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and for dynamic equilibrium 


(m]fA<n + [d](Aq] + [k] fAq] = 0 


(3-19) 


where, for example 


[m] » [A] T [I 1 ] S [A] . (3.20) 

All matrices are defined in Appendix A. 

Solution of the aeroelastic problem, Eqs. (3*l8) and (3.19), is 
accomplished in two stages. First, the nonlinear equations for the 
static equilibrium position, Eq. ( 3 - 18 ) , are solved using a Newton- 
Raphson scheme, described in Appendix E. This method has proven it- 
self as one of the best solution techniques available in geometric- 
ally nonlinear analyses [31]. It is extremely accurate and possesses 
second-order convergence. Furthermore, the completely numerical 
formulation of the FEM, together with the small number of modal co- 
ordinates, reduces the cost of forming and inverting the Jacobian 
matrix, Eq. (E.3) in the appendix. This cost usually would have to 
be accounted for as one of the major drawbacks of the Newton-Raphson 
method. Another drawback is its sensitivity to the choice of the 
initial solution vector. 

In the present case, the nonlinear solution from the previous 
value of pitch setting, 8, is taken as the initial guess. For the 
first value of 8 the linear solution is chosen as the initial 
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guess. The iteration is. considered to have converged when the abso- 
lute change of each generalized coordinate during an iteration cycle 
is less than 10 4 . If the method fails to converge, the system 
usually can be considered as being statically unstable [ 57 ]. Such a 
situation was not encountered for the flap-lag problem in hover. 

Next, the dynamic equations, Eq. (3.19), are formed and converted 
into first-order state variable form, 



- Cm] 1 [a] 1 - [m]"“ [k] 


-l 


Aq 


HI I [ 0 ] 


% 


(3.21) 


Assuming solutions of the form 



( 3 . 22 ) 


results in a standard eigenvalue problem 

(A]{y) * A(y] ( 3 . 23 ) 

where [A] is a constant nonsymmetric matrix. Equation (3.23) is 
easily solved using one of the available eigenvalue solvers. The 
eigenvalues appear in complex conjugate pairs 


\ - k t l“ k • 


( 5 . 24 ) 


Thus, the perturbed motion about the static equilibrium position is 
stable when all ^ < 0. The stability boundary is obtained by 
systematically varying the pitch setting 0 until one of the ^ is 
zero. More details can be found in Section 5.1. 
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SECTION 4 


APPLICATION OP THE METHOD 
TO THE FLAP- LAG- TORS ION AEROELASTIC PROBLEM 
IN FORWARD FLIGHT 


In this section the Galerkin finite element method (GFEM) 
is applied to a more practical problem, namely, the coupled flap- lag- 
torsion aeroelastic stability and response calculation of hingeless 
rotor blades in forward flight . Hie torsional degree of freedom 
generates an additional equation of motion, associated with this added 
degree of freedom. The torsional degree of freedom also yields a 
number of nonlinear bending-torsion coupling terms in the structural 
operator. Since the torsional equation of motion is of second order 
with respect to the spatial variable, as compared to fourth order for 
bending, special care has to be exercised when discretizing the 
equations via the Galerkin finite element method. When solving the 
discretized dynamic equations, a major complication arises due to the 
forward flight condition which introduces periodic time varying 
coefficients in the equations of motion. Due to this flight condition 
a large number of additional aerodynamic loading terms appear in the 
equations of motion. 

4.1 Brief Description of the Equations of Motion 


The equations of motion for the flap-lag- tors ion problem in 


forward flight ara coupled noollnear, nonconeervatlve, partial differ- 
ential equationa with periodic coefficients, tte structural operator 

is taken from Reference (561. tte Inertia end aerodynamic loads we 
taiten from Reference [2 ] . 

The geometry of the problem la described in Figures 3 and 

^ additlon to the ^flunptions made in Section 3.1 for the hover 
problem, the following assumptions are made : 

1. The blade has an angle of built-in twist Q (x ) 

BO' 

occurring about the undeformed elastic axis. Recall 
that the undeformed elastic axis is assumed to be 
straight and coincident with the feathering axis. 

2. Die blade croas-3ectional aerodynamic center, center 
cf gravity, and elastic center are distinct points. 

The tension center coincides with the elastic center, 
i.e., x jj * 0. 

3. The elastic torsional deformations of the blade occur 
about the deformed elastic axis. Root torsional 
deformation, due to pitch link or control system flex- 
ibility, is assumed to occur about the feathering 
axis. 

1*. Crou-Mctionu rtiffnene and inert!. prop.rt.las, off- 
sets, and airfoil chord vary along the blade. 

5 • Structural or mechanical damping of viscous type is 
included. 
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6. Aerodynamic effects associated with forward flight 
introduce cyclic pitch variations, thus, the pitch 
angle is given by 

# r (*) - « 0 + 8^ »ln * * 8 lo co. * , 
and the total geometric pitch angle is 

°G ( *0 * * °B (i 0 ) + °r (l|r) * 

7. The inflow is represented by the following general 
functional form 

X(x , if) ” Xq k^(x) + ^ (x) sin t 

+ X Q kj(x) cos * . (4.3) 

8. The effect of reversed flow is included in an exact 
manner. 

It should be pointed out that the influence of axial forces 
on the torsional rigidity of the rotor blade and the effect of cross- 
sectional warping due to torsion is neglected in Reference [ 56 ]. 
Furthermore, the effects of stall and compressibility are not included 
ir the aerodynamic loads of Reference [2 ] . Although the above effects 
may be important for certain rotor blades and certain flight condi- 
tions, no attempt is made to include them in the present study, since 
its primary objective is the application of the Galerkin finite 


(4.1) 


(4.2 ) 


element method 


In applying the ordering scheme, according to Assumption 5 
in Section 3.1, the following orders of magnitude are assigned to the 
various parameters in this study: 




i 22 

R ' / 


li 

' / 


- v’ 0(1) 


H , sin t , cos y , — — , * 0(1) , 


®0 ' ®ls ' «lc - 0( V > 


I . t> . V , x. , \ - o(e_) 


I « I 

l * 


“o ' 2 ’ “o ^ 


o(e-) . 


Using these assumptions, the coupled equations of motion 
for forward flight are presented below. Note that the equations are 
written in the reference frame e, which represents the undeformed 
blade. 


Axial equilibrium: 


T ,x * Pxl 


Lag equilibrium : 


- (M, + w - v *) , 


q 3I»x + P yl + P yA + P yD 0 


Flap equilibrium : 


(M, - 5JA v + w T) 

2,x Y ,x ,xx ,x ,x 


+ q 2I,x + P Z I + W P zD = ° ‘ 


Torsion equilibrium : 


M z,x + \ ^11 + ScA + ScD 


The corresponding boundary conditions are : 


at x Q = 0: 


v = w = v =»w =* 0 

• X ,x 


at x Q = 1: 


- \x ' ®*„ ",xx * 5 ,x f ' «3I,x * 0 


(4.5) 


(4.6) 

(4.7) 


(4.8a) 

(4.8b) 

(4.9a) 
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+ Gj(j> V +wf+Q - 

d > x ,x fXX ,x x ,2 I,x 


0 


(4.9b) 


M 3 * “ ^2 * M x = 0 

5 « o 

■me boundary conditions at the free end, ^ - l , are natural bound- 
ary conditions, expressing the fact that the shears, moments, and 
tension at the blade tip are zero. At the blade root, J * 0, the 
boundary conditions for bending involve only geometric quantities, 
i.e., the bending displacements and slopes. The mixed boundary condi- 
tion for torsion. Equation (4.8b), is a result of the root torsional 
spring. 

Equations (4.4) - (4.9) are written in a general form which is 
most suitable when using the Galerkin finite element method to dis- 
cretize the spatial dependence. All quantities appearing in these 
equations are defined below. 

The elastic moments are given by: 


(4.9c) 

(4.9d) 



(SI, cos 2 R. 8 g + J! } sin 2 R o 8 g )t ^ + (St, „SI 3 ) 

[( 2 w ,xx'* v ,n )sln2H 0 V**,xx cos 2E C , (4.10a) 



- (SI, -Elj)[(i f „ + .In 2R, «, + 

- (21, sin 2 R c 0 G * SI 3 cos 2 R o # G » ^ 



cos 2R Oj 
c G 

, (4.10b) 


48 


r 1 /-2 

2 ^ v ,xx 


- w ) sin 2R 
,xx c 


V w 
,xx ,xx 


cos 2R 

c 


« G ] , (4.10c) 


M 

X 




. (4.10d) 


The distributed force and moment vectors, per unit length of 
the undeformed elastic axis, are expressed as: 


E = 

AAA 

pe+pe+pe , 
x x y y z z 

(4.11a) 

a 33 

AAA 

ae+q_e+q_e 
be x T y ^z z 

(4.11b) 


In general, these loads contain inertia, aerodynamic, and structural 
damp ing contributions, denoted by the subscripts I, A and D, 
respectively. In writing the equations of motion, (4.4) - (4.7), the 
final form of the loads, Reference [2], has been used. Note, however, 
that due to a more consistent application of the ordering scheme, the 
torsional inertia load, Equation (4.13c), differs in some higher-order 
terms . 

Inertia loads : 


p xl = ^*1 + *0 + 2 ^ ' 

•4 • • •• 

p = m(v - v + 2 (5 w - 2u + x T cos 8_ + x T ©_ sin ©_) , 

yl p 1 u 1 u u 


(4.12a) 

(4.12b) 
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(4.12c) 


p zl * " V e i + X Q ) - W - 2p p V - x z © G cos © G ) 


l 3l " ’ “ *i cos ® G (i l + 5 0 } + d G (i , 


G m2 


W Sin 2 ®g 


- Si 


< 3 2 I 3 S ^ Sin ° 0 (5 1 + V - ^ Si - 2 ® 0 + V 3 C ° s2 ® G ‘ 


= Si 


q ll * m ijtsin © G (v - v + v ^ + £ Q )) 

. cos « 0 ( - (0 p ♦ 5^ ♦ i 0 ) - 5 . 2Pp j 

+ 4( ' ' f 1 + + 5 o> ' 2 V ?» - (1^ + 1 B} ) 

• % * * * » , V - - 1 BJ ) 


• [cos 29 G (i - s p r x * 2j, T x ) + sin 29 n ( i + j J] 

2 - 2 


,x‘ 


- ^ « G + I m3 cos‘ 9 g ) 25 x (i + ; jjt ) 


- T .* Vi + ”,X Vi + Vi 


Aerodynamic loads : 


P yA = ~^ (& G F l" F 2 + ®G V F 2 + ~T 4 + P p ^® G F 1 “ 2F 2 )▼ 

n 

+ (8 G V + 2 -f F 1> F 5 \* + ‘ e G P 1 - 2F 2 + ® G h )Fj 5 x 

' F 1 F 2 * + (8 0 F 1 * 2F 2 - 2e p Fj)» ” + (9- )v 


.(4.13a) 


.(4.13b) 


(4.13c) 


w 

.X 
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- E? w w +(PF_v + F„F, v + F.F- w )A 
3 ,x ,x pi 2 3 ,x 13 jx'* 

+ ( -2v + F. v )F, w w + (f. v + F? v )w A 

1 ,x 3 ,x ,x 1 3 ,x' ,x * 

C 

+ (8 0 F 2 + 2 -T F l* - (# G F 1 - 2F 2 F 2 F U ♦ 

+ (# G F 3 V + P 2 + * + (# G F 3 \x - 2P p 8 * 2F 3 "x * F 1 ** 

• • 

+ F Z F, w A + F, <i> w v + ( -2v ¥ +F-V w + F, A v )w 

3 4 ,x 3 >x v 1 ^ T 


• * - * ? y V i 

,x X ,x ,x 3 * x ‘ 


+ Q^wv-ww + w A + | i v] , 


(4.i4a) 


■zk 


r » - ®G F ! + F 2 - ®G V F X + P p F 1 J + (F 2 - 28 G F 1 )F 3 \x 

+ F 1 F 3 \x ’ F 1 * * F 1 J \x + <*3 ' 'l*. *,x 

- 2F 1 F 3 J ,x * + < F 2 ' 28 g V" + F 1 * -' F A * 

• • • • 

+ (F, w - 2F. <j»)v + F, v W + V W] , (4.l4h) 

J }X 1 J jX 


• * 


ScA 


: - r b 2 (0*5 - x A )(l - x A )(© G + <i»)F 1 

- F * 1 ^2 - e G F 1 )F I + 6 P F 1 7 + < F 2 - 28 G F 1 )F 3 J ,x 

+ f i f 3 *,x - 4 * + f i 8 “ X + (I 3 - 4* X \x 

- 2F X F 3 J ,x * + (F 2 - 28 G F l* + F 1 * 

• * * • • 

F, w v - 2F A v + F, v w + v w] , (4.14c) 

J fX 1 J fX 


where 
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(4.i4d) 


F 1 * ®i + *0 + 14 1 sin * 


x 1 + 4 7 P P cos * 

(4.l4e) 

R 

u — cos y , 

(4.l4f ) 

S 1 (1*5 - i A ) 

(4.i4g) 

-2 R 2 

b j (0.5 - x A )(l - i A ) . 

(4.i4h) 


Damping loads : 


P yD 

“ %L V ' 

(4.15a) 

P zD 

• 

‘ SSF^ ' 

(4.15b) 

ScD = 

• 

S S t ^ • 

(4.15c) 


The coupled flap- lag-tors ion problem is thus defined by Equa- 
tions (4.5) ” (4.7), together with the elastic moments. Equations 
(4.10), and the loads, Equations (4.12) - (4.15). Note, that the 
tension T will be el imina ted using the axial equation and correspond- 
ing boundary condition. Equations (4.4) and (4.9d). The axial dis- 
placement, u, will be replaced using the assumption that the blade 
is inextenaional in the axial direction, an assumption which is commonly 
made in rotary-wing aeroelasticity. 
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^.2 Discretization of the Equations of Motion Using the 
Gale r kin Finite Element Method 

The first step in solving the equations of motion, presented 
in the previous section, is the discretization of the spatial depend- 
ence. This is accomplished through application of the Galerkin finite 
element method. Subsequently, modal analysis is used to reduce the 
number of discrete unknowns describing the problem. 

The procedure followed here is similar to that described in 
Section 3.2, for the hover problem. Therefore, only the major steps 
will be outlined. However, special emphasis i3 placed on the appropri- 
ate modeling of the torsional degree of freedom. 

The approximate global solution given by Equation (3.7) is 
extended to include the torsional deformation: 

V • 

3 xM Mx 1 

This solution is now substituted into the flap-lag-torsion equations 
of motion. Equations (^.5) - (^*7), and the corresponding boundary 
conditions. Recall, that in the extended Galerkin method the shape 
functions $ need to satisfy only the geometric boundary conditions. 
Therefore, both the natural boundary conditions at the blade tip, 
Equations (^.9a-c), and the mixed boundary condition, due to the root- 
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torsional spring Equation (*.8b), contribute to the boundary residual 

The weighted Galerhln residual, obtained through appropriate 
combination of the weighted differential equation and boundary condi- 
tion residuals is given below. 


fl 

J [• ] 

J Q m 



„ + Gu^> w - v T) - a 

_ + v + v J) + n 

*' X T ,X fXX W ,x lJ ,X + *2 


M 


Xj x 


P T +p. + -n v 6 

yl p yA p yB 


J SI * PtA * PeD 

V ’ll* <W SrD 



* ‘W* 


(M 3 ,r + 5 J *,r i ',xx-\x T+ « 3 lJ -»3 


' ( V + V* 5 , x ** 


- M 



+ [ Voi' T 


(^•17) 


M x _ 


x 0 =° 


5 ^ 



Integrating by parts and cancelling the boundary terms, the 
final expression, corresponding to Equation (3.10), iu obtained as: 




0*.l8) 


In the interior of the element the displacements are given by: 
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- C?(i e )](a e (y)} 


(V.19) 


For bending, the same cubic Hermite interpolation polynomials 
as in the hover cane are used, see Appendix C - The nodal parameters 
are the lag and flap displacements and slopes at the element boundar- 
ies, see Figure 3b . This element satisfies the requirement of 
continuity of the global solution, since it provides interelement 
continuity for bending displacements and slopes. Bending strains vary 
linearly within the element which goes beyond the mini mu m requirement 
of constant strain within the element. 

The torsion equation of motion is of second order with respect 
to the spatial variable. Thus, a linear interpolation will achieve 
the required C Q continuity and constant strain. However, in the 
coupled bending-torsion analysis, it is desirable to use a torsion 
element which provides the same accuracy as the bending element. This 
allows discretization of the torsional variable with the same number 
of elements as needed for the adequate modeling of bending. 

In the present analysis, an improved torsion element, providing 
linear variation of torsional strain, is obtained by using the torsion 
deformation at the element mid-point as additional nodal parameter, see 
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Figure 5. Thus, N - 3, and the are quadratic interpolation 

polynomials ( 58 ], given in Appendix C. 

In general, refined finite elements can be obtained using any 
number of internal nodes. An alternative approach is the use of higher 
order derivatives (second for bending and first for torsion, or higher) 
as nodal parameters. These higher-order elements, however, experience 
difficulties in modeling concentrated loads. Furthermore, the bound- 
ary conditions involving the higher-order derivatives must be satis- 
fied. Therefore, such elements were not considered. Results for free 
vibrations of beams, using elements with internal nodes or higher- 
order continuity can be found in References [59] and [60]. 

In conclusion, it can be stated that the elements selected in 
the present study are the most basic (or simple) elements which yield 
a consistent for mula tion for coupled bending and torsion. This takes on 
an additional significance in light of the large number of nonlinear 
terms which have to be modeled. The exact form of the element inter- 
polation polynomials 7, J] , and ^ is given in Appendix C. 

The element displacements, Equation (^.19), are now extended 
over the global domain in the sense of Equation (2.7) and then sub- 
stituted into the integrated Galerkin residual, given by Equation 
(U. 18 ). This yields the nonlinear, periodic element equations. 

([ij] + [I^(a e )]){ a e ] 

+ ( [Cj] + (D|] + [D*l + [C* (a* )] + [T^(a e )] - 
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- - («;(»')] * ic^t.')] . [d' 2 (4')] 

+ ( 1 tJ ?V» - 1 Hi 1 ) 

i-efl i-1 

+ ([b'I . d'] t (k^I . [aJ] + £aj(«®)] + (ig(. e )] 

t (*|(a e )] + tAj(a')] j [a')t (F*J + {F*} + IirHa 1 ) 8^ 

" £ for e»l,2,...,E . (4.20) 

All element matrices in Equation (4.20), indicated by the 
superscript e, are defined in Appendix F. The structural operator 
is associated with the matrices [bJ] and [B*j. The axial tension 
results in the contributions represented by [t£], [T^] and 

The inertia loads are included in [ij], [l£], [c*], [cj] , [r®], 

( F j ].* ^ C Ax^ 411(1 ^ C Ax^' where the last two matrices are due 
to the axial shortening effect. The aerodynamic loads are contained 
lath, [d'], [d|], (d' 2 ], (d*], [D^J, [aJ], [aJI, U'] and 

tP'l matrices. The structural damping effect is represented by 
[Dg]. Finally, the [B c ] matrix accounts for the root torsional 
condition, where the Kronecker delta, 8 gl , indicates that this term 
is only present in the first element, i.e., the element at the root of 
the blade. 
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The functional dependence of the element matrices on the nodal 
displacements is as indicated in Equation (4.20). Note, that the 
matrices in Equation (4.20 ) have both single and double numerical sub- 
scripts. The first subscript is an identifier of nonlinear terms. A 
first subscript having a value of 2 or 3 is indicative of quad- 
ratic or cubic terms, respectively. A second subscript is attached to 
all velocity-dependent element matrices. All element matrices are 
evaluated using six-point Gaussian quadrature. The nonuniform element 
properties are included in the numerical integration. 

Next, the element matrices are assembled into the complete 
system matrices. The nodal parameters within the nonlinear element 
matrices are replaced by their modal representation, 

(a) * [A] [q] , (4.21) 

using 

shapes of the rotating blade. Subsequently, the modal reduction 
process is completed by pre- and post-multiplying the system matrices 
with the modal transformation matrix, [A], and its transpose. For 
more details regarding the treatment of nonlinear terms, see Appendix 
E. 

The final equations of motion, in terms of the reduced set of 
M modal degrees of freedom, can he written symbolically as : 

G = [m(q)][ q } + [d(q, qJH q } + [k(q)]{ q } + ( f } * 0 . 

(4.22) 


lag, Mj, flap, and torsion free vibration mode 
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All matrices in Equation (4.22) are defined in Appendix G. Note, that 
the inertia, damping, and 3tiffness terms have both linear and nonlinear 
contributions. Also recall, that for forward flight, most matrices 
have periodic coefficients, i.e.. 


£(q,q,$,t) *= G(q , q , \ , ♦ + 2t) , 


with the common period being 2 t, which corresponds to one blade 
revolution. 

For convenient numerical treatment. Equation (4.22) i 3 re- 
written in first-order state variable form. 


{y} 


( -Mq)] 1 ((d(q, <j)]{q} + [k(q)](q] + (f) ) 

! / , (4.2?) 

( ) 


and 


(y) 




(4.23a) 


It should be understood that 


(y) ~ ['}(%,*)} 


So 


is a periodic function depending nonline arly on {y} . 


4.3 Method of Solution 

As indicated before. Equations ( 4 . 22 ), governing the aero- 
elastic behavior of an isolated blade, are a set of coupled, nonlinear 
ord in a r y differential equations with periodic coefficients. It has 
been shown in previous research, References [l] and [6l], that in 
forward flight the aeroelastic stability and response problem of an 
isolated blade is strongly coupled with the overall equilibrium of 
the helicopter. The overall equilibrium of the helicopter in forward 
flight is normally obtained by enforcing the overall force and moment 
equilibrium of the complete helicopter. Such an analysis is called 
trim analysis. The results of this trim analysis are used as input 
to the aeroelastic analysis of the blade. Details of the trim proced- 
ures and the mutual interaction between the trim analysis and the 
aeroelastic analysis are presented in Section 4.3.3. 

A solution to Equation (4.22) will provide both stability 
and response information. A response solution of Equation (4.23), 
subject to known initial conditions, can be obtained by numerical 
integration. However, this approach can become expensive in terms of 
computer time, since the transient response might be significant for 
a large number of rotor revolutions. Furthermore, numerical integra- 
tion has proved itself to be a somewhat unreliable tool when dealing 
with periodic systems, unless it is complemented by explicit stability 
information. Hierefore, in this study, the desired information about 



th ' “ erMll “ tlC Stablllt i' « f «• M— ** obtained ft. the eigendata 
extracted from a linearized system. The linearized equations are 
obtained by linearizing the nonlinear equations. Equation (t. 23 ), 
about a suitable equilibrium position. An equilibrium position 
considered to be suitable if perturbations about it are sufficiently 
small, such that terms containing nonlinear products of the perturbs- 
tlon quantities can be neglected. For the case of forward flight, 
the primary source of excitation is periodic, Urns, the approximate 
nonlinear steady-state response of the blade is used as an equilibrium 
position about which the equations are linearized. Subsequently, the 

stability analysis of the linearized perturbation equations is explic- 
itly considered. 


4 * 3 * 1 ^5 linea r Steady-State Reside of ftrirvH,. 

Systems 

Previous analyses, see References [2] and ( 61 ], have used a 

simple harmonic balance technique to calculate the time dependent 

equilibrium position, ttis approach is algebraically very tedious. 

Furthermore, since terms above the first harmonic ,-er- not ret ained , 

the accuracy of th. resulting equilibrium position is somewhat 

questionable. Another shortcoming of the harmonic balance method is 

its failure to provide sufficient information on the stability of the 
response solution. 

In the present study, the equilibrium position is obtained by 
applying an algorithm which has been developed in Reference [62] for 
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the steady-state response calculation of linear systems. In order to 
apply this method to the aeroelastic problem in forward flight, 
quasilinearization [63], [64] is used to transform the nonlinear equa- 
tions (4.23) into a sequence of linear equations. The steady-3tate 
response for the linearized system at each iteration step is then 
obtained using the method of Reference [62]. In the limit, the 
sequence of linear problems converges to the solution of the original 
nonlinear problem. Additional and general information on quasi- 
linearization can be found in Reference [64] . 

The method employed in this study has two key ingredients, 
namely, quasilinearization and the linear periodic response solution. 
These ingredients are described in detail below. 

Quas ilinear izat ion , as formulated in Reference [63 J, employs 
the first-order equations, i.e., Equation (4.23). However, when the 
Galerkin finite element method is used for spatial discretization it 
is more convenient to start with the second-order system. Equation 
(4.22). These equations are expanded in a Taylor series about a 
previous solution, keeping only linear terms. The linearized second- 
order system is then transformed into first-order form. 

til 

During the K iteration step, the previous solution is 
K— X 

denoted by [q (*)}. This solution has to satisfy the boundary 
conditions exactly, i.e., it must be periodic, 

(q K " 1 (t)] » (q K ” 1 (y + 27T] , 
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and its derivatives with respect to t must be known up to the 
highest order appearing in the derivatives of G (see Equations 
(^.25) below). Furthermore, it is assumed that the derivatives of 
G exist and that they axe also periodic in Thus, 

& 3g 3g 

2 + ^ «a + + ~ 8*a + ••• 

^ h $ 

(4.24) 

To write Equation (4.2b) in a more compact fashion the following 
quantities are defined; for additional details, see Appendix G. 






3 + {/} 


(4.25a) 


{F K } = Cd(q, K ~ 1 , a K ’ 1 )lC4 K “ 1 } 


+ 

ik(a K * 1 )H? K ' 1 } * tfj , 

(4.25b) 

v 3 g 

(s 1 * * 


(4.25c) 

[c x ] = -5 

(0 s - 1 , f-1) 

(4.2 5d) 

vr & 

[jfi = — 
% 

(a 11 - 1 ) - imia 11 - 1 ;] . 

(4.2 5e) 
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The variation, 5^, which is the difference between the current 
A.nri the previous solution, i.e., 

* (<i K i - ii"" 1 ) . 

is now substituted, together with Equations (4.25), into Equation 
(U.24), which is rewritten as: 

[ifift*} + [C K Uq K ] + [S K Uq K ] 

- [C K ]C4 K-1 } " [S K ][q K-1 } + (F*} - 0 . (4.27) 


Equation (4.27) represents the linearized second-order system during 
the K th quasilinearization step. Using 



and (y*" 1 } 



9 


Equation (4.27) is transformed into first-order from: 


(y K l = [aV 1 ,^,*)]^) f tb K ( X K_1 , ♦)} , (U. 28 ) 


where 


(A K ] 


.[mYMc*] 1 


[I] 


•lift' 1 [s K ] 


[0] 


( 4 . 28 a) 
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The variation, 8^, which is the difference between the current 
and the previous solution, i.e., 


% - Cq K ] - [q K_1 ] , 


(4.26) 


is now substituted, together with Equations (4.25), into Equation 
(4.24), which is rewritten as: 

(M K ] {q K } + [C K ]{a K ] + [S K ](q K } 

- [C K ](q K ” 1 } - [S K ]{q K_1 } + (/} = 0 . (4.27) 

Equation (4.27) represents the linearized second-order system during 
the quasilinearization step. Using 




.K 

a 


K 

a 


. f K-l, 

and [y } 


.K-l 

a 


K-l 


a 


Equation (4.27) is transformed into first-order from: 


{/} = [A K ( } r K “ 1 , ^ K ' 1 , t)Uy K } *• (b &, (jr K ” 1 , t)} , (4.28) 


where 


[A K ] 


.[/r 1 [c K i ; 


[I] 


■lift 


[S K ] 


Lu] 


(4.28a ) 



(U.28b) 


/ -[M K ]" 1 (-[C K lt4 K * 1 } - [S K ](q K - 1 } ♦ f/l) 

(b K l « 

( {0 } 

Both [A K ] and (b K 1 are periodic in ♦ (with period 2n) and 

K-l 

dependent on the previous solution, £ 

The solution of Equation (U.28), i.e., the periodic response 
(y (♦)}, is then used as the previous solution for the next iteration 
step. This process is continued until convergence is achieved, at 
which step (y (*)) represents the periodic steady-state response of 
the nonlinear system, which is denoted as (y ( ♦ )} . This solution is 

used as the equilibrium position for the stability calculat ion . 

% 

For periodic systems of the type encountered in rotary-wing 
aeroelasticity, quasilinearization is particularly attractive because 
the numerical treatment required for the linear analysis of such 
systems is available. Quasilinearization is a generalized Kewton- 
Raphson method. It possesses second-order convergence [6M . Selec- 
tion of the initial solution is crucial to the success of the method. 
Additional information on quasilinearization and its properties, as 
applied to rotary-wing problems, can be found in Reference [63]. 

In the present study, two options for initiating the quasi- 
linearization process were implemented. Stability information in 
forward flight is usually plotted as a function of the advance ratio 
4. Therefore, at a given value of 4, either the linear response 
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of the system (with all nonlinear terms deleted) was used as the 
Initial solution, or the nonlinear response for a previous lower 
value of 4 was used as the initial solution. 

The linear, periodic response solution of Equation (4.28) is 
calculated using the algorithm developed in Reference [62] for linear 
systems. Recall that Equation (4.28) is the linear approxima- 

tion to the original nonlinear problem, Equation (4.23); therefore, 
the method of [62] is applicable. Hie solution procedure is outlined 
below, using the notation associated with the quasilinearization 

step. 

Urabe [65] has shown that if the multipliers, , k = 
1,2,...,2M, of the homogeneous part of Equation (4.28), 

(y K ] - [A K (i)](y K ) (4.29) 


are all different from one, then Equation (4.28) has one and only one 
periodic solution; of period 2*, given by 


(/(*» 


[» K (*)] 




)] 


• l (b K (.)lds * (y K (0)l! 


(<*.30) 


with the initial condition 

r 2 IT 

(y K (0)} - ([I] - [^(SX)] r 1 [® K (2ir)] / [ < ! ,K (s)][b K (s))ds, 

J 0 

(4.31) 
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where [$ (♦)] is the fundamental matrix of Equation (4.29), defined 
by 

[**(*)] * CA K (y)][^ K (^)] and [^(0)] - [I] . (4.J2) 

K 

It ia important to note that (y (O)} exists only if ([I] - 
(® K (2 t)]) ia nonsingular. This is the case when the magnitudes of 
all the characteristic multipliers of Equation (4.29) are different 

1. If all |Aj < 1 , then the homogenous 
system, Equation (4.29), is asymptotically stable, according to Fluquet 
theory [66], and (y K (*)} from Equation (4.30) is tne desired steady- 
state solution, lhis response is then used as the previous solution 
in the next iteration step. If any one |.AjJ > 1, then the homogen- 
eous system is asymptotically unstable. In this case, although 
Equation (4.30) Is still mathematically valid, the periodic solution is 
physically meaningless. 

To find the periodic solution, first the initial conditions, 
Equation (4.3l), are calculated. The transition matrix at the end of 
one period is evaluated using Hsu's approximate, semi analytical method 
(see Reference [67]). The periodic matrix [A (t)] is approximated 
by a series of step functions. For each step, the matrix exponential 
is approximated by a finite number of terms of its defining series. 

The integral In Equation (4.31) is evaluated by taking a constant 
value of the Integrand in each step and performing an ordinary summa- 
tion. Thus, both [ ■*> (2 tt ) ] and (y (0)) are evaluated simultaneously 


from one, i.e., |A^| f 
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during a single pass over one period. 

K-l 

Provided that the system linearized about {y (t)} is 
stable, the linearized system. Equation (^.28), is integrated using 
(y (O/j, i.e., Equation (^.Jl), as the initial condition. A fourth- 

order Runge-Kutta integration scheme with Gill coefficients and con- 
stant stepsize is used. This implies that Equation (U.30) is not used 
at all in the computations . 

Instead of storing the response [y (*)] over one period (for 
use in the next quasilinearization step), a Fourier analysis with a 
finite number of terms is performed. The periodicity of the response 
is checked by integrating over several periods (i.e., blade revolu- 
tions ) until the Fourier coefficients obtained in two subsequent 
periods agree within a desired accuracy. 

The flow chart presented in Figure o summarizes the steady- 
state response calculation of periodic, nonlinear systems, as formu- 
lated in this section. The index K used in Figure 6 denotes the 
quasilinearization iteration index which is set to zero or one, depend- 
ing on whether a linear or nonlinear initial solution is used. The 
integer NCONV indicates convergence of the response solution 
(y (t)), as compared to the solution from the previcxis quasilineari- 
zation step. At the beginning, NCONV is set to zero. When the 
response has converged it is set to one (see Step 10 in Figure 6). The 
quantity £ is a small prescribed number used in the convergence test 
and in the periodicity check. The integer K is the maximum 
allowed number of quasilinearization iteration steps. Note that 
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Figure 6 also includes the stability determination of the final linear- 
ized system, which is described in the next section. 

4.3.2 Stability of the Linearized Periodic System 

Stability is determined by deriving linearized perturbation 
equations about the equilibrium position, (y(t)}. Computation of this 
time-dependent, nonlinear equilibrium position has been described in 
the previous section; thus, at the present stage, the vector (y(t)} 
is assumed to be known. Again, when using the finite element method, 
it is more convenient to deal initially with the second-order system. 
Therefore, let 

[q] * (q} + {£q} . (^-33) 


The equations of motion. Equation (4.22), are now expanded in a Taylor 
series about [q] . Since, by definition, the perturbations (£q] are 
small, only linear terror are retained. Thus, 


8G 

a# 


G + 


& 

+ 





(4.34) 


Recognizing that G(£, ([,£,♦) * 0, since (j is the solution of 
Equation (4.22), Equation (4.34) is rewritten in first-order form as: 
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where 


(4.35b) 


IS] 

[C] 

[M] 


3g . _ 

is > 

^0 


& 

~ ( a* l) 
*1 



(4.36a) 

(4.36b) 

(4.36c) 


The ' tabUlty 0f the ‘ bov ' bOBOg.n.ou., Umar, periodic ayatem 
1. dctcmload according to Floquct theory [66], Th. tradition oatrla 
of Equation (4.35b) U defined by 


[*(<')) = [A(*)][$(*)] ^ r*( 0 )j = 


(4.37) 
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The characteristic multipliers are denoted by 

\ " \r + 1 * \i » k - 1,2, ...,2M , (4.38) 

where M is the total number of modes used in the analysis. Knowing 
the , the characteristic exponents 

\ " ^ + 1 w k 0**59) 

are determined from 

k ‘ 2T57T 10 ( 4b* 4i> > O' -WO 

1 -1 / \t v 

“k " ~j tan v 7 ^ ) * (4.40b) 

The linearized aeroelaatic system is asymptotically stable if all 
l\l < 1 > 1 * e *> 5jj < 0 • If any one |A^| > 1 , i.e., ^ > 0, 

the system is asymptotically unstable. Stability boundaries, repre- 
senting neutral stability, are obtained when |AjJ - l or ^ - 0. 

Comparing Equations (4.3C) with Equations (4.25c -e), it is 
evident that the periodic matrix [A(£ , v , *)] in the linearized 
stability equations. Equation (4.551) ), is identical to [A X ] in the 
response equations, Equation (4.28), when ^ K ’ 1 is replaced by £. 

This means that after a converged response solution £ is obtained, 
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calculation of the transition matrix and the characteristic exponents, 
i.e., Steps k and 5 in Figure 6, is repeated one more time in order to 
determine the stability of the linearized system governed by Equation 
(^•55b). The transition matrix at the end of one period is evaluated 
using the aemianalytical numerical scheme presented in Reference [67] . 

k. 3.3 Trim Procedures 

It has been shown in References [l] and [6l] that, due to the 
inherently nonlinear nature of rot ary -wing aeroelasticity, the correct 
treatment of the forward flight aeroelastic problem requires a coupled 
treatment of the dynamics of the blade and the overall equilibrium of 
the entire helicopter. In practice, overall equilibrium of a h* j.i- 
copter in forward flight is accomplished by pilot manipulation of the - 
controls, such that the control variables yield a trimmed flight con- 
dition. The mathematical procedure used to simulate the equilibrium 
of a helicopter in forward flight is commonly denoted by the term 
trim analysis. 

A simplified trim analysis, suitable for aeroelastic applica- 
tions, has been developed in Reference [6l] . The present study is 
baaed upon somewhat improved trim procedures, which are presented in 
Reference (63) • These trim procedures are suitable primarily for 
helicopters employing hingeleas rotors. For the sake of completeness, 
the main features and assumptions of the trim procedures, described in 
detail in Reference [63], are summarized below. 
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1. Blade elastic flap dynamics, based upon a single mode 
representation axe included. 

2. Blade precone, P , offset x A , and linear built-in 
twist are accounted for. 

3. Variable inflow, given by Equation (**.3)» and reversed 
flow effects are included. 

k. The helicopter is assumed to be in straight, steady flight 
at constant speed, i.e., u - 0. 

5- Rotor shaft dynamics, i / 0 , and tail rotor effects 
are not considered. 

6. Stall and compressibility effects are neglected. 

7. The linear steady-state flap response is evaluated 
using the method described in Section U.3.1. 

Two different trim procedures are considered. 

Propulsive Trim (see Figure 7 ), simulates actual forward 
flight conditions. The advance ratio u and the weight coefficient 
(approximately equal to the thrust coefficient C^) are 
specifiec and horizontal and vertical force equilibrium, as well as 
r?ll and pitch moment equilibrium, are satisfied. The computed trim 
outputs are the rotor angle of attack, Op , the thrust coefficient 
C^, , the pitch components 0 Q , , C^ c , the inflow X , and the 

steady-state flap response. 


Moment. Trim (see Figure 7 )» simulates conditions under which 

a rotor would be tested in a wind tunnel. Pitching and rolling 

moments are maintained at zero. Force equilibrium is satisfied 

implicitly, since the rotor is mounted on a supporting structure. 

When using this trim procedure, 4, 0 Q , and 0^ are specified and 

C . 0 .0 , and the flap response are computed. 

T Is lc 

The control settings, ©q, ©^ s , and the inflow A 

obtained from the trim analysis are then used as input in the aero- 
elastic (single blade) analysis . 

The iterative process, combining trim and aeroelastic analysis 
(see also. Figure 8), can now be described as follows: For a specific 

blade a given flight condition (Steps 1 and 3 in Figure 8 ), 

a) a trim analysis is performed to compute the trim values 

6 9 , 9 and (During the first trim analysis - 

O' Is' lc 

S^p U in Figure 8 - a simplified elastic flap response is 
calculated within trim. ) 

b) Using zhe trim values from a), the aeroelastic analysis is 
performed to calculate the stability and steady-state 
response of the rotor blade (Step 5 in Figure 8 ). 

c) The flap response from a) is compared with the more accur- 
ate flap response obtained in b ) . If they agree within a 
desired accuracy, the results obtained in b ) are the solu- 
tion of the aeroelastic problem and the computations are 


terminated. 


V 


d) If the convergence test in c) is not met, the flap response 
from b) is imposed on the trim analysis. Subsequently, 
the trim values are recomputed (Step 7 in Figure 8). These 
new, improved trim ve *s are then used to reevaluate the 
steady-state flap resp^. ^ j from the aeroelaatic analysis. 
This process is continued, iteratively, until the converg- 
ence test indicated in c) is passed. 
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SECTION 5 


RESULTS AND DISCUSSION 


The numerical results presented in this section are intended 
to illustrate the application of the Galerkin finite element method to 
rotary-wing aeroelasticity. It should be noted that the main emphasis 
is on the application of the method to bingeless rotored helicopters. 
However, the method is equally applicable to other types of rotors, 
such as teetering or articulated rotors. Furthermore, the method is 
eminently suitable for analysis of modern bearingless flexbeam-type 
rotors which have a complicated redundant structure. In the following 
section, three separate grcxips of results will be presented: First, 

results for some rotating beam-type free vibration problems will be 
given in terms of free vibration frequencies and mode shapes. Next, 
the coupled flap- lag aeroelastic stability boundaries in hover will be 
evaluated. In this case, the critical collective pitch angle at which 
instability occurs will be presented as a function of the inplane 
frequency. Finally, results illustrating the aeroelastic stability 
and response of the coupled flap-lag dynamics of a hingeless rotor 
blade in forward flight are presented. For this case the governing 
parameter is the advance ratio and the results will illustrate the 
variations of system damping and response as a function of advance 
ratio. 


77 



ORIGIN* l •••* '* 

OF POOR OUA'J 

5 .1 Some Computational Details 

Due to the different nature of the hover and the forward flight 
cases, two separate computer programs were written. This approach 
permitted a more efficient treatment of each problem. 

In the hover case all matrices are time invariant. Thus, for 
a uniform blade the basic parameters of the problem, i.e., the lag and 
flap stiffness, the collective pitch angle and the inflow, could be 
extracted from the element matrices. Hence, these quantities appear 
explicitly as scalar factors in the finite element equations. There- 
fore, after selecting the appropriate number of elements, the element 
matrices need to be calculated and assembled into the system matrices 
only once. Next, for a particular value of lag and flap stiffness, 
the uncoupled free vibration modes are calculated and the modal reduc- 
tion is carried out. The reduced linear and nonlinear system matrices 
are stored. For a particular value of pitch setting the reduced system 
matrices are multiplied by the basic parameters and added up to give 
the final equations, Equations (3.18) and (3.19)* The iterative 
Newt on-Raphs on solution to the nonlinear static equilibrium equations, 
Equation (3.18), is considered to have converged when the absolute 
change of each generalized coordinate is less than 10 

The search for the stability boundary is started from zero 
pitch setting. The pitch setting is then increased in equal steps 
until the system becomes unstable. To obtain a more accurate estimate 
of ® c , bisection is used in the interval in which the system has 
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become unstable. Subsequently, the pitch setting is incremented 
further to determine whether the system remains unstable. 2his search 
algorithm routine is a modified version of that used in Reference [57]. 

In forward flight most matrices contain periodic coefficients. 
Thus, calculation of the element matrices, assembly of the system 
matrices, and modal reduction has to be repeated at each value of the 
dimensionless time ♦ . Efficient programming and use of efficient 
methods to compute stability and response is, therefore, crucial for 
the effective treatment of this problem. 


The efficient numerical computation of transition matrices is 
very important. Transition matrices are evaluated using the semi- 
analytical scheme presented in Reference [67]. To find the transition 
matrix at the end of one revolution (i.c., common period) and to 
compute the integral in the initial conditions, the interval 0 to 2 t 
is divided into N ^ equal sectors. In each sector the required 
matrix exponential is approximated by the first five terms of its 
defining series expansion. 


th 

The periodic response of the K linearized system is comput- 
ed with a fourth-order Runge-Kutta scheme with Gill coefficients. The 
first-order system is integrated over several revolutions using a 

step size of 2 t/n . During each revolution, N , where N = 
iKi rev rev 

1 denotes the interval 0 to 2tt, the response is Fourier analyzed, 
keeping harmonics, i,e., (l + 2N^) terms, for each generalized 

coordinate. An absolute change of less than 10”^ for each Fourier 
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coefficient indicates the periodicity and the convergence of the non- 
linear response; see Steps 8 and 9 in Figure 8. The Runge-Kutta 
integration and Fourier analysis routines are improved versions of 
those used in Reference [62 ] . 

The trim results computed in this study are generated using 
the trim program from Reference [63]. In order to perform the trim 
and aeroelastic analyses iteratively, the trim program had to be 
linked to the finite element aeroelastic analysis program. 

The eigenvalues of the dynamic system in hover. Equation 
( 3 . 23 ), and the transition matrices at the end of one period (forward 
flight case) are obtained by utilizing the IMSL Math emetics! Sub- 
routine Library provided by IBM. First, the matrix is preconditioned 
by reducing its norm through exact diagonal similarity transformations. 
The matrix is then reduced to upper Hessenberg form by orthogonal 
similarity transformations. FinSLlly, all eigenvalues are computed by 
using the QR method with origin shifts at each iteration. 

All element matrices are evaluated using six-point Gaussian 
quadrature . 

Results based on the global Gale r kin method are computed using 
a modified version of Power's computer code. Reference [57]. Five 
nonrotating modes for each flap, lag, and torsion are used to obtain 
the uncoupled free vibration frequencies and mode shapes of the 
rotating blade. 
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5 .2 Free Vibration Results 


In order to validate the finite element method, uncoupled free 
vibrations of a cantilevered beam were considered first. 

The convergence properties of the selected finite elements can 
best be evaluated by computing the free vibration bending and torsion 
frequencies of a nonrotating uniform beam. Figure 9 shows the relative 
accuracy of the finite element solutions, as compared to the exact 
solutions, when the number of elements is increased. The first bending 
and torsion frequencies, curves labeled IB and IT respectively, are 
obtained with excellent accuracy even when only two elements are used. 
The second and third bending and torsion frequencies (curves labeled 
2B, 2T and 3B, 3T) are obtained with less than one percent error 
when three and five elements are used, respectively. Furthermore, it 
is apparent that the cubic interpolation bending element end the 
queuiratic interpolation torsion element provide approximately the same 
accuracy. For comparison, the results for the first torsional fre- 
quency, when using the linear interpolation torsion element, are also 
shown. Die performance of this element is considerably inferior when 
compared to the refined torsion element based upon quadratic interpo- 
lation. All elements exhibit uniform convergence. This was antic- 
ipated, since the finite element model for this conservative problem 
can be derived from a variational principle. 

Results for the first bending frequency of a nonrotating, 
nonuniform beam are shown in Figure ]0. The beam has linear width and 
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depth taper, the width and depth values at the free end being 80 
percent of the respective root values. The ensuing quadratic mass and 
cubic bending stiffness distributions are approximated by using the 
functional values at ten equidistant points and assuming linear varia- 
tion in between. Based on this approximation, the mass and stiffness 
properties within the elements sure modeled in three different ways. 
First, the mass and stiffness within the element are assumed to be 
constant, using the average of the respective nodal values. Second, 
mass and stiffness are as aimed to vary linearly within the element. 
Third, the beam properties are integrated together with the element 
shape functions, using six-point Gaussian quadrature. All frequencies 
are referenced to a six-element solution, where the exact functional 
form of the beam properties is included in the element integrals. 

Figure 10 shews that numerical integration of beam properties yields 
accuracy comparable to the case of a uniform beam (see curve IB in 
Figure 9 ). Modeling the beam properties as linear within the element 
gives acceptable results. However, it is computationally more tedious. 
Lastly, using the constant model would require a significantly larger 
number of elements, solely for the purpose of modeling the nonuniform 
mass and stiffness properties. 

In the present study, numerical integration of the nonuniform 
beam properties is chosen. It gives the best results and is computa- 
tionally moat easily implemented, since the element matrices are eval- 
uated using numerical integration. 

The free vibration problem of a uniform, rotating beam, having 
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In- and cxit-of- plane bending and torsion degrees of freedom, is 
considered next. Hie fundamental uncoupled lag, flap, and torsion 
frequencies are calculated using tvo methods: 

1. A global Galerkin method in which five nonrotating modes 
of a uniform beam sure used to obtain the fundamental, 
uncoupled, flap, lag, and torsion frequencies of a 
rotating beam, and 

2. The Galerkin finite element method of weighted residuals. 

Comparison between the two sets of results, showing the un- 
coupled rotating flap and lag frequencies, w and <2 , and the 

r -L LI 

torsion frequency, as a function of the nonrotating frequencies 

is shown in Figures 11 end 12, respectively. The results for the 
fundamental frequencies are identical when three or six elements are 
used, indicating that good accuracy can be obtained with a relatively 
small number of elements. 

A comparison of the first- and second-flap mode shape, obtained 
using both methods is presented in Figure 13 for the nonrotating case 
and two different speeds of rotation. For this case, eight elements 
were used because it was more convenient to plot the results when 
slightly more elements were used. From a numerical convergence point 
of view, a smaller number of elements wcxild have been adequate. 

From these results it is clear that excellent agreement between 
the global Galerkin and the Galerkin finite element method is obtained 
for both frequencies and mode shapes. 
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The free vibration results presented In this section show that 
the selected cubic interpolation bending and quadratic interpolation 
torsion elements sre very accurate, even when only a small number of 
elements is used. It is reasonable to expect that these elements will 
also perform very well In the ae roe las tic analysis, since the aero- 
dynamic loads do not depend on the strains. Modeling of strains, i.e., 
the higher-order derivatives of the elastic degrees of freedom, 
usually governs the accuracy and convergence of the finite element 
solution. 

The actual number of elements used in the aeroelastic analysis 
will be governed by the number of modes retained in the modal reduction 
process . 

5.3 Re stilts for the Case of Hover 

In this section the Galerkin finite element method is applied 
to a typical rotary-wing aeroelastic problem. The coupled flap-lag 
aeroelastic problem in hover [l], [68] is a convenient and simple 
example which can be used to illustrate the Galerkin finite element 
method. The finite element equations for this problem have been 
derived in Section 3*1. 

In calculating numerical results, certain simplifying assump- 
tions were made, since the objective of this section is primarily to 
illustrate the application of the Galerkin- type finite element method 
to rotary-wing aeroelasticity . 

These simplifying assumptions are listed below: 
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1. The inflow vu assumed to be constant over the blade and 
equal to the inflow at 0.75 span, i.e., 

X - (ao/l6)[(l ♦ 2U0/a a) 1 ^ 2 - l] . (5.l) 

2 . Hub and tip losses were not included . 

5. Structural damping was assumed to be zero. 

Pertinent values of the data used in the calculations are 
presented in Table 5*1 below. 

TABLE 5.1 

CONFIGURATION PARAMETERS FOR FLAP-LAG IN HOVER 



The numerical accuracy of the method can be best seen by 
comparing a global Gale r kin solution, baaed on one uncoupled 
rotating elastic mode for each degree of freedom, with a local 
Galerkin finite element solution in which the blade is 
represented by three finite elements and where, for consistency 
with the global method, one uncoupled elastic rotating mode is used 
for each degree of freedom to reduce the number of nodal degrees of 
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freedom. The comparison between the two methods is presented in 
Table 1. For two separate collective pitch angles ©, all pertinent 
values associated with these cases were evaluated. The agreement 
between the two methods is excellent when considering that only three 
elements are used to represent the blade. Similar comparisons were 
made for a variety of other cases, the results are not presented here 
since they would liave been repetitive of the behavior illustrated by 
Table 1. The cases presented in Table 1 were stable configurations 
because the elastic coupling parameter R. whs taken as R c - l.C. 

Convergence properties of the Galerkin finite element method, 
when ipplied to the aeroelastic problem, are considered next. It is 
important to note that convergence of Galerkin-type methods in aero- 
elasticity can be established only by numerical experimentation 
In rotary-wing aeroelastic problems this is further complicated due to 
their nonlinear nature. Convergence of the method can be investigated, 
alternatively, by varying the number of elements while retaining a 
fixed number of modes in the modal reduction process, or by changing 
the number of modes and maintaining a fixed number of elements. 
Convergence with respect to the nonlinear iterative solution technique 
is not investigated here. It should be noted, however, that converg- 
ence is very rapid. At most, three iteration cycles are required. 

Figure 1 U illustrates the convergence of stability boundaries for 
three different values of the elastic ' upllng parameter when the 
number of elements is allowed to -/ary froa three to six, while the 
number of nodes used in the modal reduction process is maintained at 
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one mode for each elastic degree of freedom. In all cases shown in 
the figure, the unstable regions tend to decrease as the number of 
elements is increased. The results for E - 5 and E - b are al- 
most identical, indicating that four or five elements are sufficient 
to capture the bending dynamics of the blade. 

Figure 15 shows the convergence characteristics when the number 

of elements is maintained at four, E - 4, while the number of modes 

used in the modal analysis is allowed to vary. The curves represent 

points at which the real part of those eigenvalues associated with 

predominant lag modes, i.e., ^ and is zero. The real part 

of the eigenvalues associated with predominant flap modes, i.e., ^ 

ajyi ^ is always negative. The couplea modes of the aeroelastic 

system, Sq. ( 3 - 25 ), can be identified as predominantly lag or flap 

modes by correlating the imaginary parts of the eigenvalues with 

the frequencies of the uncoupled rotating beam vibrations, system 

ii, st 

stability boundaries are obtained by piecewise combination of 1 

Lag” ann " 2 ^ Lag" curves such that, overall, the lowest value of 

0 is maintained. For R - 0.0 md R - 0.4, system stability 
c c 

is determined by the behavior of the predominant first lag mode. The 
predominant second lag mode is always stable, i.e., s, L < 0 • For 

R - 0.0 it is interesting to note that convergence trends differ, 
c 

depending on whether the blade is -oft inplane or .tiff lnplana. 

In the Intermediate range of elastic coupling parameters, 0.5 < * c < 
0.8, the possibility of a flap-lag type Instability dominated by the 
second lag mode can occur. For R c - 0.6 the system is unstable 
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1 . 6 , the 


due to the second lag mode when w <1.6. Above = 

LlL 4JlL 

first lag mode causes the system to become unstable. For R * 0.8 

c 

the second lag mode solely determines system stability. The first lag 
mode also beccmes unstable, however, at higher values of 9. Finally, 
although not shown in the figure, it should be noted that the system 
is stable for R above 0.9 for the complete range of to.. 

C LlL 

investigated. In conclusion, it can be said that the difference be- 
tween taking two or four modes (i.e., M 3 2 V3. M 3 k) is rela- 

tively small with respect to the behavior of the first lag mode. This 
holds for values of up to 2 . 5 , which is rarely exceeded in 

practical hingeless helicopter blade configurations. However, it is 
absolutely essential to use fom modes, since the second lag mode 
governs system stability for intermediate values of elastic coupling. 

In Figure 16 the number of elements is kept constant at six, 

E 3 6 , and the number of modes is changed. The behavior of the two- 
and four-mode model follows basically the same trends as those 
observed when four elements are used (see Fig. 15 ). When using six 
modes, as compared to four moaes (i.e., M = 6 vs. M = M, the same 
results are obtained, except for the instability associated with the 
third predominant lag mode However, this instability occurs at very 
high pitch settings and is therefore not relevant for system stabil- 
ity. These results indicate that it is sufficient to use a total of 
four modes, two modes each for the elastic flap and lag degrees of 
freedom. As a cross-cneck, the four-element, four-mode solution for 

R 3 0.6 is also shown in Fig. 16. The results are essenti ally the 
c 
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same as in the six-element, four-mode solution. A difference occurs 
at fi^-0.8 and *» 1.6; however, it is very small, indicating 
that four or five elements are sufficiently accurate even when four 
modes are used. 

Figures 17 and 18 illustrate the influence of the elastic coup- 
ling parameter B c on the real parts of the eigenvalues, which are a 
measure of modal damping and determine system stability. The pitch 
setting 6 is kept constant and two different values of lag frequency 

are considered. The predominant flap modes axe very stable and 
damping associated with them remains almost constant when the elastic 
coupling is changed. Damping of the predominant lag modes changes 
considerably with elastic coupling. Damping of the first mode is 
further strongly influenced by the value of It is also evident 

that the instability of the second lag mode is relatively weak and 
usually it can be eliminated by including a small amount of structural 
damping in the analysis. 

A comparison of stability results when using uncoupled versus 
coupled mode shapes in the modal analysis is depicted in Figure 19* 

Use of coupled mode shapes should lead to more accurate results; 
however, it has the disadvantage that the coupled free vibration 
problem has to be recomputed for each increment of pitch setting. A 
total of four modes is used in each case. The four lowest frequency 
coupled modes include the two lowest lag and flap modes (predominant) 
when < l.h and, for >1.^, the lowest lag and the three 

lowest flap modes (predominant ) . This means that the actual modal 


representation changes with the inplane frequency w when simply 
using the lowest frequency coupled modes. Thus, the analysis is not 
able to capture the behavior of the second predominant lag mode. 

Figure 19 shows that accuracy is not improved when using coupled 
modes. On the other hand, computing time increased roughly four- fold 
for this example. Therefore, use of coupled modes is not recommended. 
It should be noted that when six coupled modes were used the predom- 
inant second lag mode was present for all values of end it3 be- 

havior was modeled correctly. However, computing cost increased in a 
prohibitive manner without any noticeable gain in accuracy when com- 
pared to the results using uncoupled modes. 

In order to be able to compare the GFEM with the global Galerkin 
method, Piwers ' computer code [57] was modified by excluding torsion 
and introducing the elastic coupling parameter. The GFEM program was 
set up such that it can represent the governing dynamic equations of 
motion from Reference [57] without torsion. Then both methods were 
applied to the sane problem. Figures 20 and 21 show that excellent 
agreement is obtained when using four or six uncoupled rotating modes. 
These are modeled by four or six elements in the GFEM and ten un- 
coupled nonrotating modes (five for flap and five for lag) in the 
global Galerkin method. Again, it is concluded that four elements 
and four modes will be sufficient to model simple flap-lag blade 
dynamics in hover. 

Figure 22 illustrates the effect of using different sets of non- 
linear equations of motion. In Reference [5] the global Galerkin 
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method based on four coupled rotating modes (one predominant lag, 
three predominant flap) is used. In the GFEM, four elements and four 
uncoupled rotating modes (M^ * M^, * 2 ) are used. Both, equations 
from Reference (57] (without torsion) and the present equations of 
motion are employed. Instability of the second lag mode at R c * 0.6 
is not shown. The trends established by the three solutioiu are, 
overall, the same. There are, however, marked differences in the 
actual numerical values, in particular, at R c - 0.0. This must be 
attributed to the fact that the equations are slightly different. 
Furthermore, it should be noted that the flap- lag aeroelastic problem 
in hover, in the absence of elastic coupling terms, is very sensitive 
to small, higher-order terms; thus, small differences in these terms 
can lead to noticeable differences in the results. This also indicates 
that unless the equations are derived in an identical manner, differ- 
encec can be observed for the same problem. These differences are, 
however, exaggerated by the fact that R c was taken to be zero. 

Figure 22 concludes the results presented for the flap-lag in 
hover problem. The numerical experience gained by applying the 
Galerkin finite element method to this basic rotary-wing aeroelastic 
problem provides valuable guidelines when dealing with the more prac- 
tical flap-lag in forward flight problem, for which results are pre- 
sented in the next section. 

The computing times required to generate the stability boundar- 
ies in this section were only slightly larger them those required when 
using the global Galerkin method. This might seem to be somewhat 
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surprisingj however, it has to be attributed to the fact that the 
f Jnite element program was efficiently structured to avoid unnecessary 
recomputation of element matrices, assembly of system matrices, and 
modal reduction. The exact comparison of computing times for each 
method depends, of course, on the number of elements employed in the 
finite element model. 

5.^ Results for Flap-Lag Blade Dynamics in Forward Flight 

Numerical results presented in this section deal with the coup- 
led flap-lag aeroelastic problem in forward flight. The finite element 
equations for the complete flap- lag- tors ion problem in forward flight 
have been given in Section ^ .2 . No previous finite element solutions 
for the stability and response of nonlinear, nonconservative systems 
with periodic coefficients are available. In view of the novel fea- 
tures of the present research, where a finite element solution to such 
systems is given for the first time, it was deemed appropriate to avoid, 
initially, the added complexity of dealing with the torsional degree of 
freedom. However, it should be mentioned that the torsional equation 
of motion and the additional terms associated with the torsional degree 
of freedom in the lag and flap equations, do not introduce any concept- 
ually new effects. They do, on the other hand, increase the size of 
the problem, imposing additional requirements on computer storage and 
computer time. 

The coipled flap- lag finite element equations in forward flight 
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were obtained by simply disregarding all torsional submatrices (third 
row and column in the partitioned element matrices, see Equation 
(F.l)) and deleting all torsion terms in the nonlinear bending sub- 
matrices. This simplified system vas carefully checked and found to 
be consistent with the ordering scheme. The method of solution pre- 
sented in Section ^.3 was formulated in a general manner and is di- 
rectly applicable to the flap- lag problem, as well as to the coupled 
flap- lag- torsion problem. 

In calculating numerical results two types of inflow were in- 
cluded in the analysis. 

1. Uniform inflow, where X is given by Equation (5 -2 ) and the 
cyclic inflow components, X fl and X , are set to zero; 
thus. 


M- tan 0 ^ + - 


'T 


2 N/ ^ + 4 


(5-2) 


2 . Nonuniform cyclic inflow, where X^ is given by Equation 
(5*2) and the total inflow is: 


X = p. tan Op + 



(l + 1.2x cos *) . 


(5-3) 


For calculating the inflow at p. * C.O in the moment trim procedure, 
Equation (5 - 1 ) was U6ed. 



The numerical results obtained for the flap-lag problem in for- 
ward flight are presented in two groups. 

First, results from the application of the Gale r kin finite ele- 
ment method (GFEM) are compared with previous solutions, where the 
global Galerkin method (GGM) was used. This group also contains 
results illustrating the numerical properties of the solution proced- 
ure for the discretized dynamic equations, described in Section ^.3. 
Based on the experience gained with the hover problem, described in 
Lection 5*2, three elements and, for consistency with the GGM results, 
a total of two modes were used. The elastic coupling parameter was 
set to R c = 1.0. The data for this group of results is presented in 
Table 5.2. 


TABLE 5-2 

CONFIGURATION PARAMETERS FOR FLAP-LAG IN FORWARD FLIGHT 



(Figures 

23-30) 

5 U = 1.417 ; 

(J 

FI 

= 1.007; 

E = 3 ; 

M 

= 2 ; R = 1.0 ; 

c 

b = 0.0313; 

7 

= 5*0; a = 0.05 ; 

N psi ’ N r*i ' *h 


variable . 


In all these cases, propulsive trim with a weight coefficient of 
3 0.01 and uniform inflow was used. The actual trim values, 



taken from the results of Reference [2], axe listed in Table 2. 

The second group of results deals with the convergence proper- 
ties of the Galerkin finite element method and with the influence of 
several configuration parameters on system stability and response. 
The data values used for this group of results are presented in 
Table 5*3. For the soft inplane blade, these properties axe close to 
those of the Boelkow BO- 105 hingeless rotor. 

TABLE 5-3 

CONFIGURATION PARAMETERS FOR FLAP-LAG IN FORWARD FLIGHT 
(Figures 31-39) 



The trim values were calculated using the trim procedures from Refer- 
ence [63]* In all cases, the fuselage pitching moment and the var- 
ious trim offsets were set to zero. Propulsive trim values for a 
typical value of weight coefficient, ■ 0.005, are shown in Table 
3- 

Parameters which remained unchanged for all forward flight 
results are: o A ■ 1.23 kg/m^ (0.00238 slugs/ft^); a • 2tt, C^ 0 = 
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C DP 3 °* 01; *L • °‘ 0; *U " 1 -° J «1 " P p ’ X A “ y F " °* 0j 

■ 0.0. Furthermore, blade pretwist was set to zero, © B ■ 0.0, 

and the blade properties were assumed to be uniform over the span. 

Other pertinent quantities are specified on the plots. 

Stability results, showing the comparison between the 
Gale r kin finite element method with the global Galerkin method, are 
given in Figures 2). The GGM results were taken from Reference [2], 
Figure 9^ and Reference [6l], Figure 9* In Reference [2] the same 
equations of motion as in the present study were employed. Three 
different values of torsional stiffness were considered. For the 
comparison, the results for the highest torsional stiffness, * 

15*033, were used. Reference [6l] dealt with the flap- lag problem 
in forward flight. However, the equations of motion did not include 
same of the higher-order terms retained in the present study. In 
al 1 three cases, stability results were obtained by linearizing the 
equations of motion about an approximate linear time-dependent equi- 
librium position. To simulate the harmonic balance method, as it 
was used in References [2] and [6l], only the first harmonic, (N^ * 
l) was used in generating the GFEM results. For two advance ratios, 
M. * 0.2 and p. ■ 0.U, ten harmonics were also used. To find the 
periodic response, the numerical integration was carried out over 
three blade revolutions and the Fourier coefficients, computed 
during the third revolution (N rev 3 3) were used. Hie number of 
azimuthal steps (per revolution) for both stability (N^ g ^) and 
response (N was taken as 120, which is identical to the 
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stepaize used in Reference [2], 


The real part of the characteristic exponent foi the predomi- 
nant lag degree of freedom, , versus the advance ratio, p., is 

plotted in Figure 23a. The value of is a measure of the overall 

damping associated with the first lag mode. The GFEM results exhibit 
the same trend as results from Reference [2 ] . The difference ln ^u, 
is roighly the same for all advance ratios. This difference must be 
attributed to the absence of the torsional degree of freedom in the 
GFEM results. This conclusion is further confirmed by the fact that 
the GFEM results agree very well with those from Reference [6l], up 
to an advance ratio of ^ ■ 0.2. When p. is increased above that, 
the additional higher- order nonlinear terms in the presently employed 
equations of motion become more important, so that eventually, at 
M- * 0.^, the GFIM result is closer to the result from Reference [2 ] . 

Stability results for the flap degree of freedom. Figure 23b, 
generally show the same behavior as discussed above for lag. It 
should be pointed out that results for the characteristic exponents 
from References [2] and [6l] exhibit the splitting, typical of per- 
iodic systems, for ^ from p. =* 0.2 on upward, while the GFEM 
yields split values only at an advance ratio of p. - O.U. At lower 
values of p. the characteristic exponents appear in complex conju- 
gate pairs. This difference between the two methods arises, probably 
due to the new procedure employed for solving the dynamic equations 
in the present study. 

Agreement between the Galerkin finite element method and the 
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global Galerkin method results, presented in Figures 23, should be 
Judged in light of the previously discussed differences in the blade 
model and in the procedure uned to obtain the linear time -dependent 
equilibrium position. Overall, qualitative agreement between the 
two methods ia quite good. Stability results at M- ■ 0.4, using 
four finite elements, brought the GFEM results in even better agree- 
ment with results from Reference [2], For the configuration in 
Figure 23, the blade motion was stable for all advance ratios 
considered. 

When comparing GFEM results based on one harmonic (N^ - l) 
with those using ten harmonics (N^ ■ 10 ), it is abvious that at 
p. * 0.2 results for ^ (Figure 2 3a) and ^ (Figure 23b) are 
almost identical. At p * 0.4 there is a minor difference for the 
real part of the lag characteristic exponent. The flap degree of 
freedom, on the other hand, exhibits a rather remarkable change in 
when using ten instead of Just one harmonic. This last result 
indicates that it is important to use more than one harmonic when 
considering advance ratios above p = 0.2. 

Figures 24a and 24b illustrate the effect of the number of 
harmonics used in the Fourier analysis on the first approximation 
(K * l) to the nonlinear steady-state lag and flap response . The 
blade tip displacements (generalized normalized coordinates) for p * 
0.4, corresponding to the stability results in Figures 23, are 
plotted during the third blade revolution. Again, the contribution of 
the higher harmonics is more signii'icant in flap than in lag. This, 



as a matter of fact, explains the remarkable change In flap damping, 
as seen in Figure 25b. 

Table 4 contains the Fourier coefficients of the response as 
plotted in Figures 24. From the relative magnitude of these coeffi- 
cients it is apparent that, while it is essential to retain more than 
the first harmonics, it probably would be sufficient to use the first 
four harmonics. This conclusion was also found to be valid for other 
configurations . 

Since the CRJ time required for the Fourier analysis was very 
small, ten harmonics (}« ■ 10 ) were used in al 1 s’ibsequent 

u 

calculations . 

Figures 25a and 25b address the question of how many blade 
revolutions, i.e., periods, should be integrated in order to obtain a 
periodic response. These plot6 show that the response is almost the 
same during the first (N rey “ l) and the second (N rey " 2) blade 
revolution. The difference between the response from the second and 
third blade revolution cannot be distinguished on the plots. 

If the maximum absolute change of any Fcxirier coefficient, as 
compared to its value from the previous revolution, is taken as a 
measure of convergence, both the second and the third revolution 
yield Fowler coefficients which have converged within an accuracy of 
10 ^ . Thus, on error control parameter based on such an accuracy 
(10 seems to be adequate. Also recall that in the derivation of 
the equations of motion the displacements were assumed to be of order 


99 


Ep m 0.2 and terms of O(e^) were neglected, as compared to 
terms of 0(l). Therefore, this error control quantity Is also log- 
ically consistent with the ordering scheme. Furthermore, It should 
be noted that the Initial conditions used for the numerical integra- 
tion theoretically insure a periodic response. The effect of approx- 
imations and numerical errors in the actual calculation of the initial 
conditions will most likely be corrected with the integration over the 
second or third blade revolution. Any further integration will merely 
lead to an accumulation of errors in the integration routine. 

Stability results (not shown), based on linearization about the 
steady-state response from the second and third revolution, did not 
change in any significant manner. This further confirms that an 
accuracy of 10 ^ for the absolute change of the Fourier ^coefficients 
is sufficient. 

Table 5 shows the first-order state variable vector at the 
first three full blade revolutions, for the response plotted in Fig- 
ures 25- The elements of this vector are the lag and flap blade tip 
speeds and displacements in the rotating system. At the azimuthal 
angle * ■ 0, this vector is identical to the initial conditions. 

From ♦ - 2tt on, only small changes occur. A monotonic convergence 
trend from one complete revolution to the next cannot be observed. 
Hcwever, by using the error c -'.trol parameter associated with the 
harmonic components, as discussed previously in this section, the 
periodicity can be accurately determined. 

Convergence of the results indicative of stability. Figures 
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26, and of response, Figures 27, with the number of quasi linearization 
steps was studied nert. As expected, results for stability and 
response change markedly when going from the linear stability and 
response solution, K ■ 0, to the first nonlinear solution, K - 1. 
Stability results baaed on the first nonlinear equilibrium position, 
J.e., curves labeled K ■ 2 in Figures 26, show a further change only 
at the higher end of the advance ratio range, i.e., at y ■ 0.-*. This 
result, again, emphasizes the importance of retaining higher order 
nonlinear terms when advance ratios above y ■ 0.2 are considered. 
Therefore, stability information should be obtained by linearization 
about a nonlinear time-dependent equilibrium position. For the 
configuration considered in Figures 26 (dat« from Table 5»2), the 
nonlinear terms decrease the stability margin of the Lag degree of 
freedom for all values of y. For the flap degree of freedom, stabil- 
ity is decreased only at y ■ 0.4; at lower values of y the non- 
linear terms are stabilizing. It is interesting to note that flap 
stability results are most sensitive to the number of quasllineariza- 
tion steps used in the analysis, i.e., nonlinearities are more 
significant for flap. As a matter of fact, at y ■ 0.4 the question 
arises if, indeed, the lower of the two valuer of the characteristic 
exponent for flap has converged with the second iteration step; see 
Figure 26b. 
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Response results in Figures. 27 show that nonlinear terms (K * 
l) increase the amplitudes of both the lag and flap response, as 
compared to the linear solution (K = 0). When performing an addi- 
tional iteration step, i.e., going from K * 1 to K = 2, the 
response remains practically unchanged. The corresponding response 
curves cannot be distinguished within the accuracy of Figures 27. 
Inspecting the Fcxirier coefficients at K = 2 for convergence, as 
compared with the previous iteration step K * 1, 1* was seen that 

the maximum absolute change of any Fourier coefficient at P ■ O.k 

—3 

is less than 10 . 

The above question of convergence of the characteristic ex- 
ponents with the number of quAsilinearization steps was pursued 
further by allowing a third quasilinearization step (K = 3)» for 
stability only, i.e., by linearizing the system about the nonlinear 
response from the second iteration step (K = 2). Figure 28 shows 
the relative change of the characteristic exponents versus the number 
of quasilinearization steps for the most critical cases presented in 
Figures 26. At P = O.b- the exact solution was assumed to have been 
obtained with the third iteration (K * 3), while at p = 0.2, the 
second iterative solution (K a 2 ) was used as reference. The 
results in Figure 28 Indicate that at p = 0.4 two iteration steps 
are needed, while at the lower value of advance ratio, P = 0.2, one 
iteration 3tep is sufficient to obtain converged values for the char- 
acteristic exponents. Thus, it can be concluded that at higher 
advance ratios, use of a nonlinear equilibrium position i3 necessary 
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to obtain accurate stability information. 

Fran the results presented, in Figures 26 through 28, the 
following procedure was established to obtain converged nonlinear 
response and stability results. First, the quasilinearization iter- 
ations are carried on until the maximum absolute change of any 
Fourier coefficient of the response is less than 10~\ Second, the 
system is linearized about this converged nonlinear steady-state 
response to yield the final explicit stability information, i.e., the 
real part of the characteristic exponents. It shcxild be pointed out 
that convergence of the characteristic exponents with the number of 
iteration steps can also be monitored, since this information is 
available at each quasilinearization step. However, computationally, 
it would be very difficult to implement a c divergence test based on 
the characteristic exponents, because their identification is compli- 
cated by the fact that the imaginary parts are not uniquely known and 
that they do not always appear in complex conjugate pairs. The iden- 
tification would be particularly difficult in the case where more 
than one mode per elastic degree of freedom is considered. 


The effect of the number of azimuthal stations in calculating 
stability and the initial conditions, N 811(1 number of steps 

per blade revolution in the Runge-Kutta integration, N ric jL» is 
illustrated in Figures 29 and 50. Note that all results presented 
for forward flight are obtained using the same step size for the 


stability and response calculation, i.e., 
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The relative change of the real part of the characteristic 
exponents at h = O.U versus the number of azimuthal steps, is 
shown in Figure 29- The solution with N^ s ^ = 120 was assumed to 
be correct. Note that the flap degree of freedom does not have 
complex conjugate characteristic exponents, i.e., two distinct real 
parts are associated with it. Overall, results besed on fourty steps 
are in excellent agreement. Even twenty steps, N p Si 3 20, give 
acceptable results. These conclusions, however, should be viewed as 
dependent upon the particular configuration considered and the number 
of modes used to reduce the nodal degrees of freedom. When ten steps 
were used, the numerical procedure broke down. Note also, that for 
the particular configuration considered here (data from Table 5*2), 
the flap results converge slower than the lag results. The reason 
for this apparently being that higher harmonics are more significant 
for the flap response than for the lag response; see Figures 2k. The 
higher harmonic response contributions obviously require a smaller 
stepaize for a certain, desired accuracy. 

Figures 30a and 30b compare the lag and flap response, respect- 
ively, for 20 and 120 steps (per revolution). Agreement between 
the response curves for the different number of azimuthal steps is 
quite good, although the absolute values of the Fourier coefficients 
differs by more than 10 ^ . When using steps, the response could 

not, be distinguished from the curves based on 120 steps, within the 
accuracy of the plot (result not shown). For this case (N^ gi * i*0), 
the absolute difference for the Fourier coefficients of the response 
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from the two different step sizes was less than 10~^, indicating 
that forty steps would be sufficient. 

Figure 30b concludes the results intended to illustrate the 
effects of the various parameters, associated with the solution pro- 
cedure, on the blade dynamics in forward flight. These parametric 
studies were necessary since previous numerical experience was limited 
and restricted to the case where the linear steady-state response was 
used as the equilibrium position. No solutions to nonlinear aero- 
e las tic periodic systems using quasilinearization were available. 

From the numerical experience gained in this study, the fol- 
lowing conclusions were drawn. It is important to retain more than 
the first harmonic in the Fourier analysis of the response. In sub- 
sequent calculations, ten harmonics, N^ = 10, were used. The non- 
linear response was considered to be periodic and to have converged 
when the maximum absolute change of each Fourier coefficient was less 
than 10~^. According to this criterion, periodicity was achieved 
with the second or third blade revolution and the quasilinearization 
procedure converged with the first (K * l) or second (K * 2 ) non- 
linear response, depending on the value of the advance ratio. Forty 
azimuthal steps were sufficient to obtain accurate response and stab- 
ility solutions. In all subsequent calculations, sixty steps, N pgi * 
N rk ^ * 60, were used. This should be adequate even when more ele- 
ments and mode shapes are used in the analysis. 

The configuration for which these results were calculated 
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(data In Tables 5*2 and 3) correspond to a relatively high loading 
case (Cy * O.Ol). The choice of parameters made above should thus 
be considered conservative when cases with a more moderate blade 
loading are sought. Therefore, the parametric study above was not 
repeated for the other blade configurations and discretization models 
considered in this study. 

The convergence properties of the Gale r kin finite element 
method are considered next. This is accomplished by changing the 
number of elements or the number of mode shapes in the modal reduc- 
tion process. All the results are based on the configuration para- 
meters given in Table 5.3. 

The relative change in the real part of the characteristic 
exponent versus the number of elements is shown in Figure 31 for the 
soft in-plane blade, = 0.732, and elastic coupling R, * 0.6. 

The number of modes was kept constant at two. As reference, the five- 
element solution was used. It is apparent that excellent convergence 
is achieved, in particular, when considering that the results in 
Figure 31 are for a high advance ratio, p. = O.k. Interestingly, the 
accuracy for the flap degree of freedom is much higher than that for 
lag. This must be attributed to the lower stability margin for lag; 
see Figures 32 . Overall, the three-element solution can be considered 
sufficiently accurate. It should be kept in mind, however, that the 
configuration in Figure 31 is stable. For a more critical case, more 
elements might be required to model the system accurately. Finally, 
it is interesting to compare Figure 31 with the accuracy for the 
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first bending frequency of a nonrotating beam in Figure 9» Aa ex- 
pected, the solution of the aeroelastic problem does require a larger 
number of elements than the free vibration problem. 

Figures 32 show system stability when changing the number of 
modes from two to four, while keeping the number of elements constant 
at E * The aeroelastic damping for the fundamental modes, 

1L 

and remains unchanged when using four modes as compared to two 

modes. The damping, i.e., real part of the characteristic exponents, 
for both predominant flap modes is practically constant for all ad- 
vance ratios. The absolute value of is somewhat lower than that 

for the first flap mode, however, both modes are strongly 

damped. The first lag mode has its lowest damping values at moderate 
advance ratios, p = 0.1 -0.2. When the advance ratio is increased, 
the results indicate that the forward flight aerodynamics have a 
stabilizing effect. Overall, the smallest stability margin occurs at 
hover, p = 0.0, for the second predominant lag mode. However, with 
increasing advance ratios, more aeroelastic damping is fed into the 
second lag mode. For advance ratios, p > 0.2, the real parts of 

and £g L , are 

roughly the same. 

The lag and flap blade tip response was plotted in Figures 33 
and 3 1 *. The configuration considered is the same as that for which 
stability results were presented in Figures 31 and 32 . 

Figures 33 show the response for three different advance ration 


the characteristic exponents for both lag modes, ^ 
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obtained by using four elements and two modes. The time dependence 
of the response in forward flight (m- ■ 0.2, OA) as compared to 
the static response in hover (m- * 0.0) is clearly illustrated. 

The time dependent contribution to the lag displacements is basic- 
ally a one per rev motion, while the major contribution to the flap 
displacements is a two per rev motion, l.e., it stems from the sec- 
ond harmonic. This is the same behavior as encountered for the 
stiff in- plane blade in Figures 2k. Another interesting aspect of 
the response curves in Figures 33 is that the displacements at the 
advance ratio p. * 0.2 are lower than those for hover. From the 
trim data in Table 3 it is obvious that there is a direct relation 
between the value of collective pitch setting ©^ and the magnitude 
of the response. The largest displacements occur at the advance 
ratio p. * OA which has the largest value of 0^. 

The effect of the number of modes retained in the modal 
reduction process on the blade response at p. * OA is illustrated 
in Figures }k. The response of the second lag mode, hg* is very 
small. Only in the reversed flow region can it be distinguished 
from zero. The second flap mode response, g^, is more Important . 
The behavior of the response associated with the first lag and flap 
mode, respectively, varies accordingly. For lag, the response of the 
first mode does not change significantly when going from two to four 
modes in the analysis. For flap, on the other hand, there is a size- 
able change. When the response of the first and second flap mode 
(from the four-mode analysis, M = M are added together, its 
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maximum value is roughly 8 percent larger than that of the responae 
based on the two-mode analysis (M *2). 

The effect of the number of modes used in the analysis was 
further investigated by considering stability of a stiff in-plane 
blade, u> u - 1.417, with elastic coupling R c =» 0.8. Results for 
the real part of the characteristic exponents are obtained by using 
two and four modes. In both cases, the blade is represented by four 
elements. The stability curves in Figures 35 exhibit the same gener- 
al behavior as encountered for the soft in-plane blade (Figures 32). 
There are, however, two important differences. The second lag mode 
is unstable at M- - 0.0, i.e., is positive. Thereafter, the 

forward flight aerodynamics introduce a considerable amount of damp- 
ing, so that at M- = 0.1 the second lag mode is more stable th an 
the first lag mode by a factor of five. A further increase in the 
advance ratio, changes the value of ^ such that it approaches 
the value of £g L , and at li« 0.4 they are practically the same. 

The other interesting point is that at n - 0.4 only the four-mode 
solution exhibits splitting of the characteristic exponents (real 
part) associated with the first flap mode. The two-mode solution 
does not capture this effect. 

Results presented in Figures 34 and 35 indicate that for both 
response and stability it is important to retain four modes in the 
analysis. Recall, that in the hover case the second predominant lag 
mode itself was the cause for system instability at certain values of 
the lag frequency <3,^; see Figure 15. In the forward flight case 
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such an Instability was not observed; however, the presence of the 
second lag and flap mode did change the response and, in certain 
cases, the stability behavior associated with the predominant first 
flap mode. Although this change did not result in a critical condi- 
tion, it certainly will have an effect on blade bending moments and 
vibration levels. This, in turn, effects the fatigue life of the 
blades and the vibration levels in the fuselage. 

The effect of the elastic coupling parameter R on the stab- 

c 

ility of the first ai second predominant lag mode is shown in 
Figures 36 and 37, respectively, for the soft in-plane blade. The 
stability margin of the first lag mode (Figure 36) increases propor- 
tionally with the value of R c throughout the entire range of advance 
- ratios. The least stable configuration is obtained for zero elastic 
coupling, at low advance ratios (u - 0.1 -0.2). The behavior of the 
second lag mode, (Figure 37), is quite different. The variation of 
damping versus the advance ratio depends strongly on the value of 
V ^t^* tt should be mentioned that the predominant flap modes 
are very stable and the damping associated with them remained almost 
constant when the elastic coupling was changed. The same observation 
was made for the hover case; see Figures 17 and 18. 

A comparison of stability results for two values of the weight 
coefficient is presented in Figures 38. The damping associated with 
the second lag and flap mode is practically unaffected by the value of 
C W r0r the first fla P mode the higher value of reduces damping 
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increasingly with the advance ratio. Also note that at U * 0.1* the 
blade subjected to higher loads (C^ ■ O.Ol) exhibits splitting of 
the real parts of the characteristic exponents for Damping 

values for the first lag mode show a reversed trend. The lower value 
of C ^-0 .005 decreases ^ considerably, in particular, at high 
advance ratios. Additional results (not presented here) for a still 
lower value of the weight coefficient, C,^ * 0 .0025, showed that the 
small stability margin of the predominant first lag mode at 1 ■ 0.1 
was reduced even further. This is consistent with the physical be- 
havior of the flap- lag instability in hover, since lower values of 
at low advance ratios result in reduction of the aerodynamic 
damping available in the lag degree of free dan. 

The results presented so far were based on the uniform inflow 
model, given by Equation (5.2). Figure 39 illustrates that the in- 
fluence of the nonuniform inflow, Equation (5*3), on the stability 
results is minor. A small difference for the lag stability results 
can be observed only at M- * 0.1. The flap stability results (not 
shown) were unaffected at all. These results suggest that either the 
configuration considered in Figure 39 is not sensitive to the inflow 
model, or the inflow representation must be more sophisticated than 
Equation (5*3) in order to have any effect on system stability. 

Finally, it should be pointed out that all forward flight 
results presented in this study were obtained by using the linear 
equilibirum position as an initial guess in the quasilinearization 
procedure. For a test case, system response and stability were 
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computed, Incrementing p. in four equal steps from p ■ 0.0 to 
P ■ 0.^. The linear response at the current value of p, as well as 
the converged nonlinear response from the previous lower value of p, 
were used as the initial guess for the nonlinear response solution. 
Both, the nonlinear response and the stability values of the linear- 
ized system, converged to the same solution for all values of p, 
regardless of which type of initial guess was employed. Since using 
the linear initial guess was computationally cheaper, this option was 
chosen. It alio has the additional advantage that results can be 
computed for any arbitrary value of advance ratio, without having to 
sweep the entire range of pi's, starting from zero up to the desired 
value of the advance ratio. 

The results presented in this section for the flap-lag problem 
in forward flight indicate that this is basically a stable configura- 
tion. Figure 35a illustrates that the second lag mode instability 
encountered in the hover problem (Section 5«2) does not persist in 
the forward flight region. Fran results in Reference [2] it becomes 
clear that inclusion of the torisonal degree of freedom introduces 
instability, usually associated with the lag motion. However, only 
for torsional ly soft blades does this instability manifest itself in 
the range of advance ratios considered in the present study. 

The computational times required for the flap- lag problem in 
forward flight are quite significant. They depend on a number of 
parameters, such as the number of ex-oents and modes, number of 
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azimuthal atepa, number of quasilinearization iterations, and number 
of revolutiona integrated. Furthermore, the configuration parameters 
also play a role, in as much aa they determine the degree of nonlin- 
earity of the system. To find the converged, nonlinear periodic 
response and linearized stability for one value of advance ratio, 
approximately 30 CRJ seconds were needed in the three-element two- 
mode case. When using four elements, hO CRJ seconds were required. 
In the case of four elements and four modes, this value increased to 
approximately 100 CRJ seconds. The computation of the element 
matrices of aerodynamic origin, which must be performed for each value 
of azimuth angle, takes up roughly 50 percent of these CRJ times. 

A direct comparison with computing times necessary for the global 
Galerkin method was not possible since those results were generated 
using the harmonic balance method to compute the equilibrium position. 
All forward flight results were generated on an 134 3033 computer. 

Lastly, it should be mentioned that the iterative trim - 
aeroelastic analysis, as discussed in Section ^.3*3 and indicated by 
Steps 7 and 8 in Figure 8, was not used to generate results in the 
present study. It would have increased the computing times 
considerably, without contributing to the basic objective of this 
study which is primarily aimed at the application of the finite ele- 
ment method to rotary-wing aeroelasticity. 
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SECTION 6 


CONCLUSIONS 


This study presents the formulation of a Gale r kin- type finite 
element method for nonself ad Joint, nonlinear aeroelastic rotary-wing 
problems. Fran the numerical results presented for the aeroelastic 
stability and response of hingeless helicopter rotor blades, the 
following conclusions can be drawn: 

1* Tlle Galerlcijl finite element method is a practical tool for 
formulating and solving rotary-wing aeroelastic problems. 
Since spatial discretization is applied directly to the 
partial differential equations, algebraic manipulative 
labor is reduced significantly when compared to the appli- 
cation of the global Galerkin method to similar problems. 
However, more computer time is spent in calculations. 

2. Four or five elements are sufficient to capture the bending 
dynamics of the blade with the same accuracy as the global 
Galerkin method. 

3. Normal mode transformation, combined with the Galerkin 
finite element formulation, reduces the number of nodal 
degrees of freedom significantly and enables one to deal 
efficiently with complex problems. 

For the flap- lag problem in hover it is 


essential to use 


two modes for each elastic degree of freed cn, since the 
second lag mode determines system stability for certain 
values of elastic coupling. 

5. The flap- lag problem in forward flight is basically stable. 
The lowest stability margins are associated with the lag 
degree of freedom at moderate advance ratios and low rotor 
loading. Inclusion of two modes for each elastic degree 
of freedom is necessary to determine blade response 
accurately. 

6. Nonlinear effects are important for both stability and 
response, in particular, at high advance ratios. 

7. Higher harmonic contributions to the periodic blade motion 
are significant, especially for flap stability and 
response. 

At this point it should be mentioned that a portion of the 
research presented in this dissertation has already been published 
\S9] and was presented at the Fourth European Rotorcraft and 
Powered Lift Aircraft Forum. 

Based on the experience gained frcm the flap-lag problem, the 
Gale r kin finite element formulation for the coupled flap- lag-torsion 
problem, presented in Section U, can be implemented directly. It is 
expected that four or five elements ere also adequate to model the 
coupled bending-torsional dynamics of the blada. 


115 


Finally, it should be noted that the Galerkin finite element 
method, as formulated in this study, provides a natural tool for 
dealing with the canplex structural configurations encountered in 
modern bearingleas flexbeam-type main rotors and tail rotors. 
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FIGURE 1: Rotor/Fuaelage Model 
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FIGURE 2 : Elements and Configuration ffcrameters of a 

F.lngeless Blade 
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FIGURE 6: Flow Chart Illustrating the Stability and Response 

Calculations of the Nonlinear, Periodic Aeroelastic 
System 
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FIGURE 7: Configuration for Trim Procedure-? 













HUMBER CF ELEMENTS 


FIGURE 9: Accuracy of the First, Second and Third Bending 

and Torsion Frequency Versus Nuaber of Elements 
for a Hcorotatiag Uniform Bean 
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FIGURE 10: Finite Element Modeling at a Tapered, Ncorotating 

Beam (Mtea and Stlffnese Within the Element are 
Modeled ae C exultant, Linear aad Ualng theatrical 
Integration)* Result a are for the Lowest Bending 
Frequency. 





FIGURE II : Comparison of the Methods in Calculating 

Fundamental Flap and Lag Frequencies of 
a Rotating Blade 
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FIGUKL 15 : Convergence of Flap-Lag Stability Boundaries V/hen the Number of Modes 

(a =■ 0 . 10 ; y = 5 . 0 ; w F1 = 1.15) 


j s Changed 




Corves represent points at which real part of eigenvalues associated 
with predominantly lag modes are zero. 
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FIGURE 1? Real Pert of Eigenvalue* (Eq. 3.2 3) aa a Function of 
Elastic Coupling 

(S - U, M - U, a - 0.10, 7 - 5.0, - 1.15) 
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FIGURE lfi: Peal Bart of Eigenvalues (Eq. 5.2 3) as a Function 
of Elastic Coupling 

(E • L; M ■ I*; a - 0.10; y - 5.0: 
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FIGURE 20: Comparison of Flap-Lag Stability Boundaries When Using GFEM Versus Global 
Galerkin Method (GGHF 

(o - 0.10, 7-5.0, u yi - 1.15) 

^Eqs. of Motion are those from Ref. (57) without torsion. 
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FIGURE 24b: Contribution of Higher Harmonics to 

the Flap Steady-State Response 


(N . = N . . = 120; K = 1; 
psi rki 

all other data from Table 
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rev 
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FIGURE 25b: Periodicity of Flap Steady-State Response 

for Three Consecutive Rotor Revolutions 

( Np si = N rki = 120; K = 2; all other 
data from Table 5.2) 
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ADVANCE RATIO 


FIGURE 26a: Effect of Nonlinear Terms on Lag Stability 

(N = N . = 120: N =3; all otner 
psi rki rev 

data from Table 5.2) 
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FIGURE 27b: Effect of Nonlinear Terms on Flap 

Steady-State Response 

(N = N . . * 120; N - 3; all other 
psi rki rev 

data from Table 5.2) 
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FIGURE 29: Canrergence of Real Part* of Characteristic 

Exponents as e Function of the Hunber of 
Axinuthsd Steps (* pgl “ N rki' K " 5; "rer " 
3; ^ ■ 0 A ; all other data fra* Table 5-2) 
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FIGURE 30a: Effect of the Hi 
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HUMBER UF ELEMENTS 


FIGUR1 31: Caerrergence of Flap- lag Stability When the 

Humber of Elaeents ii Changed (M ■ 2 ; 

u., - 0.732; R - 0.6; n - O.l*; propiUiTe 

Ul c 

tria; - 0 . 003 ; all other data free 
Table 5-3) 


ADVANCE RATIO 


FIGURE 32a: Lag Stability When the Number of 

Modes is Changed 

(E = 4; u = 0.732; R = 0.6; propulsive 
Li- C 

trim; C = 0.005; all other data from Table 5-3) 
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FIGURE 32b: Flap Stability When the Number of 

Modes is Changed 

(E = 4; * 0.732; R^ = 0.6; propulsive 

trim; C = 0.005; all other data from Table 5.3) 
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FIGURE 34a: Lag Steady-State Response When the 

Number of Modes is Changed 

(E = 4; U) = 0.732; R = 0.6; propulsive 
Li- C 

trim; C = 0.005; all other data from 

w 

Table 5.3) 
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ADVANCE RATIO 

FIGURE 35a: Lag Stability When the Number of 

Modes is Changed, Stiff Inplane Blade 

(E = 4; oj = 1.417; R =0.8; rronulsive 
LI c 

trim; C = 0.01; all other data from 
w 

Table 5.3) 
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ADVANCE RATIO 

FIGURE 35b: Flap Stability When the Number of 

Modes is Changed, Stiff Inplane Blade 

(E = 4; = 1.417; = 0.3; propulsive 

trim; C = 0.01; all other data from 
w 

Table 5.3) 
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FIGURE 38a: Lag Stability for Two Different 

Weight Coefficients Versus Advance Ratio 

(E = 4; M = 4; uJ _ = 1.417; R = 0.8; 

c 

propulsive trim; all other data from 
Table 5.3) 
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FIGURE 38b: 


Flap Stability for Two Different 
Weight Coefficients Versus Advance Ratio 


R = 0.8; 

— Lii C 

propulsive trim; all other data from 


(E = 4; M 


4; oj l1 = 1.417; 


Tshlo 5.3) 


ADVANCE RATIO 


FIGURE 39: Influence of Nonunifcrm Inflow on 

Lag Stability 

(E = 3, M = L, = 0.732; R = 1.0; 

LX c 

propulsive trim; C w = 0.005; all other 
data from Table 5.3) 
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TABLE 1 


COMPARISON OF GLOBAL GALERKIN AND FINITE ELEMENT GALERKIN METHODS 

a -0.05; 7 = 5-0; “ F1NR - O.h-, “iaNR “ 1#li R c = 1,0 

Finite element results based on three elements and one mode per 
degree of freedom; Global Gale r kin results based on one mode per 
degree of freedom. 



Global Galerkin 
0 -0.20 

GFEM 
0 - 0.20 

Global Galerkin 

0 - 0.45 

GFEM 
» = 0.45 

- (0 = 0) 

1.179967 

1.18026 

1.179967 

1.18026 

LI 





- (o-o) 

1.140290 

1.14150 

1.140290 

1.14150 

“fi 





<4 

0.016719 

0.016484 (-1.4$) 

0.087953 

0.086768 (-1.3$) 

0 

«i 

0.070289 

0.066521 (-5-^$) 

0.187816 

0.177959 (-5-2$) 

^1L 

-0.026198 

-0.026168 (-.11$) 

-0.066079 

-0.065838 (-.36$) 

^1F 

-0.307892 

-0.308048 (.05$) 

-0.281260 

-0.281628 (.13$) 


A$ 


GFEM - Global Galerkin 


100 


Global Galerkin 






TABLE 2 


FRORJISIVE TRIM VALUES* 



x 0 

°R 

°o 


°lc 

1 3 


.0000 

.2970 

-.0000 

.0000 

i 


.0031 

.2652 

-,0605 

.0038 

0.2 


.0187 

.2559 

-.1178 

.0072 

0.3 


.0406 

.2874 

-.1878 

.0113 

0.4 

.03939 

.0672 

.3577 

-.2853 

.01 66 


# r l 

Based on = O.Ol j taken from results of Ref. [2], Used 
for configurations where data in Table 5.2 was employed. 
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TABLE 3 


FR0HJI3IVE TRIM VALUES* 


n 

x o 

“r 

®0 

0 1B 

°lc 

0.0 

.05000 

.0000 

.1432 

-.0003 

.0002 

0.1 

.02518 

.0094 

.1084 

-.0247 

.0078 

0.2 

.02057 

.0406 

.1086 

-.0487 

.0151 

0.3 

.03463 

.0875 

.1454 

-.0867 

.0248 

0.4 

.05653 

.1250 

.2110 

-.1501 

.0382 


« 

Baaed on data in Table 5 * 3 ; C w - 0.005, uniform inflow. 
Calculated with trim program from Ref. [ 63 ]. 
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TABLE 4: mm I EH COEPT IC IENT3 0 T STEADY-STATE RESPONSE WHEN US DIO CHE AND TEN HARMONICS 
(Reaulta pertain to reaponae plotted In Flfa. 2 U ) 



NH 

- 1 


NH 

10 


h l 

8 1 


h l 

*1 


CONSTANT TERPB 


CONSTANT TERMS 

1 

-0.63047E-01 

0 . 11936 Btoo 

1 

-0.62985K-01 

0.ll936E*oo 


SINE TERMS 


SINE TERMS 

1 

0.144336*00 

-0.57833E-O1 

l 

0.144178*00 

-0.57849E-01 


COSINE 

TERMS 

2 

-0. 49644 L-oe 

0.3656OE-O2 




3 

0.88oo4e-03 

0.12577E-CE 

1 

-O.U065E-O1 

- 0 .i 5 ce 7 E .01 

u 

0.2 3353E-03 

-0.39137E-03 




5 

0.18575E-04 

-0 . 572 60E-04 



6 

-0.12173E-04 

-0.198696-04 



7 

-O.15202E-O4 

-0.10889E-04 



8 

-0.13324e-o4 

-0.96163K-05 



9 

-0. 12 404 E- 04 

- 0 . 8008 IE -05 



10 

-0.1D909E-04 

-0.719 37 * -05 




COSINE 

TERMS 



l 

-0.U636E-01 

-0.15032E-01 



2 

-0.21675E-01 

-0.30516E-01 



3 

-O.73437E-03 

-0.28532E-02 



it 

-0.22548E -l>4 

-0.222 30C -03 



5 

0.50308E-04 

-0.532 74E-04 



6 

0. 39837B-04 

-0. J3600K-04 



7 

0.41726E-04 

-0 . 3o667E-o4 



8 

0.40969E-04 

-0.29E84E-04 



9 

0.41057E-04 

-0.29574E-O4 



10 

0. 398o8e-o4 

-0.28595 E-04 


? 

O 
7 » 
CD 
: n 

c 





TABLE 5 

STEADY-STATE RESPONSE AT FIRST THREE COMPLETE BLADE REVOLUTIONS* 


* 


g x U) 


gl (t) 

0 

.15852 

-.04481 

-.09814 

.07095 

2 jt 

.13830 

-.04874 

-.09685 

.07137 

4ff 

.13822 

-.04821 

-.09740 

.07072 

6n 

.138U9 

-.04858 

-.09724 

.07loe 


* 

First-order state variable vector at * - n *2i r, 
n m 0,1,..., 3. Results pertain to response 
plotted in Figs. 25. 
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COORDINATE TRANSFORMATIONS 


Th« transformation* between the various coordinate systems are 
given below. Note that for rotations of order 0(« D ) the smau angle 
as sumption is used a 
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APPENDIX B 


INTEGRALS AND COEFFICIENTS 


y Q - T, co. 0 G - C sin © Q 
« t) sin © G + C cos © G 


li ' “ • h 

JJ ? “ ‘ *3 

p dA ■ m 

JJ p y Q dA - . xj 

If 0 z 0 “ ■ ” X I 

rr 2 

i/ p ’ « • ** 

r r 

Ji ° 52 “ ■ ^3 

p r] ^ dA * O 

J c " x o "’‘o * \ 




CO. 8 g 

.in 8 g 
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}--J * ; 0 » 2 i)di 0 

1 

T i<i 0 ) - J i(i x * ; 0 )ai c 


•r > 1 fX ° 

“ ( V 2 J o 




2 n^ b # R 

C " " ' J 

r R 


and 


2 a o, b R 

A « 


*b 


where - equivalent constant value of semichord b 



*22 " 001,2 R c °G + gI 3 

- Elg - (Eig - SIj) sin 2 R c 


R c °G 


"25 ■ k (81 2 - f ty ,in 2K c ®G 

b 33 ■ il 2 ,ln2 B C ®o * E1 3 co * 2 "c ®g 
- BIj + (Big - ETj) #ln 2 R c C G 


h ' £ < S5 2 - E1 3> “• 2E c ®g 


B o 22 ’ 2<1 rf c °* 2 *G * E b3 511,2 ®G } 
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B m23 * (1 »2 ' V 28 G 


B «33 " 2(1 »! ,ln2 *G * Vj CC> * 2 V 


V " (1 »2 ' W C °* 29 G 


B nO ’ * *m3* 


"IS ‘ ° X I sln °" 


Ijc - « Xj CO. 8 0 


Fcr uniform blade : 


r - t 

6R 


« 2 

LIRE 


-^7 < B 1 '> k 

"o n i 


-2 

“Film " J 4 ^ P 1 ^ 

“o a i 


T1KR 0 2 2. ' : 

n r(l m2 " I «3 ) 
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APPENDIX C 


ELEMENT INTERPOLATION 
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APPENDIX D 


MATRICES FOR HOVER 


The various matrix operators and element matrices required for 
the treatment of the flap- lag problem in hover are presented below. 
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^ p p ( *o + ^ 


+ r 


e( *o + V 7 x + ( 1 )2 x2 - “T *o ( *o + 2 *i ) 


e *c/*0 + 2i i } - ( *o + *1 } / x 


[F e } 


/ l*] T (F}<£ 

J Q e 


[G] 


[G 6 ] 


I- 0 


- 2P P 


2P p 


[ * (['if [o] [*])<£ 

J r\ 


Uii ■ 


[if] - / (1*] T [!][*])« 

1^0 e 


[K x ] 


» r Z e 

[K*] * J (t^] T [ ][^])dL^ e 
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[kl - [A] T ([B] s ♦ [T^ 8 - [ ^ ] S + [Ai ]S + [A 2 f)[ A ) 
[m] * [A] T [i i ] S [A] 






B„ D 4 
x 




D (T. D ) 
X 1 X 


D (T. D ) 
X 1 X 


T 1 = I [(1+ 25 1 ) - *o(*o + 


[s L ] = [A] T ([B] S + [^f - [ K] _ 1 S + [A 1 ] S )[A] 
[ 8nl ] = [A] T [A] 

[T x ] = 


T. D 
1 x 


T, D 


ir> 



I 



2v ji ( ) 0 

,x ^ 




) 





7 


~,x 



0 


0 





7 7 1 

~,x ~,x 


eO 
h 2 


/ ® 7 1 dx, 


~,x 



eO 

s 2 





o 

*** 


0 


t 
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APPENDIX E 


SOLUTION OF THE NONLINEAR EQUILIBRIUM EQUATIONS 

IN HOVER 

The final nonlinear static equilibrium position equations 
are given as 


(f) • + [s^U 0 )!^ 0 ) - (c) - o , (E.1) 

see Eq. (3*18). The solution to this system is obtained by applying 
the iterative Newt on- Raphs on technique [U7] . The solution increment 
during iteration step n can be expressed as 

■ (A, - (A,.! - <1^ <*- 2 > 


where the Jacobian matrix [j] is given by 


[J] 



!s l> * -75 • 


(E.3) 


In order to illustrate the evaluation of the nonlinear term and 
its derivative, it is helpful to use indicial notation [33]. In 
the following lower/upper case subscripts will be used for the 
element/system nodal parameters, while Greek subscripts will be used 
for the generalized modal coordinates . 
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On the 


element level the nonlinear term is, see Eo. (3-lM, 


[A^Ha 60 } 


(E.U) 


or, in indicial notation, 


e \ eO ,eO 

(a 2 >ljk • a 3 • *k 


(E-5) 


Using Eq. (3.16) to express a£° in terms of its modal representation, 


(E-5) becomes 


<#«k • *3° • v 4 


(E.6) 


^ A >«a ' ‘f ' 4 


(H.7) 


Here, (A^ A) is a new third-order tensor. This tensor is now 
assembled to yield an expression for the system, namely. 


/ 0 0 
A ^ua a J ' ^ 


(E.8) 


How, the modal reduction is completed as indicated by Eq. (3-20), i.e., 
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\l ^2 A hja A jp % *a 


(E-9) 


or 


i 


t \ 00 

s llL •>pa ^ ’ Sx 


(£. 10 ) 


Rote that expression (E.IO) can be written in matrix form which will 
yield the final expression for the nonlinear term as it is presented 
in Eq. (E.l). 

The derivative of 1s simply taken according to the 
chain rule 


3 / v 0 0 

v 0 ®RL i/k Sc 


^NlA/k 


* 1 ° 




^Nl^ijk Si + ^NlAiJ q / 


(E.11) 


“ \ 
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APPENDIX F 


ELEMENT MATRICES FOR FORWARD FLIGHT 


The element matrices for the coupled flap- lag -tors Ion problem 
in forward flight, see Equation (U.20), are defined below. Each 
coupled element matrix can be written in a partitioned form, e.g.: 




^ a ilf^ 

[a ilt ] 

1^1 ' 

^ a ifl^ 

J 

^ a iff^ 

[Ajjt] 


if 

^ a itf^ 

1 


(F.l) 


Since some of the element matrices are rather lengthy expressions, 
they will be defined in terms of the sub-matrices. Sub-matrices which 
are not listed below are Identically zero. 




* 


<*;! • r J, 


yf » p <Vi*®'2> 

1 

n^ (6 G r i- 2F 2 i 

.T 




1 

1 

♦ Vk )P 3 ! 

n 

F l*2 


t]V T p P 

T 

1 

r 

1 




— P 1 

* 5f, (p 2 - 29 G F i )f } 

1 

1 

my, r i F , I 

1 

-ai T 

•5 

di 

e 


4- 

1 

— - 

- — — - 


£ 1 ; a<£ t s p r l 

+ 4V* (F- - 2 O F )F_) 

1 

1 

1 

1 

-tf 

£ 1 ; » *1 




f * 1 ■ r J, 



(-2F, 





%*£xx 2B 2J 


[I £’ ‘ J„ * !,**:£cc ( 2 B i. * 55 *?„) 


jW^ ( -B , v® + 2B. w® ) 
— ,xx O ,xx 4 ,xx 


y n T ( 2 b, A® - gJ A® ) 

xx ~, xx 4 y ,x 7 


2,**C* 2B >5 


An T b w e 

— , xx 23 , xx 


+ A tj T Gj v* 
Sx-ix ,xx 




V CO 




a <°G F 3 J %- 2 0p ;e 
' 2t 3 * F ! 


T -f 

aa f 3 \* 


ai T f 3 f 4 5% 


dx 


ia T 5 ; a f 3 v s 


-e 

w 

*X 


+ F 1 VV tF 3* e v> 


dx 
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TJ _T 

mO x 



fK £lJ • J Q f - yf 


[K 1 TL ] " J 0 x is - fc?x (x is ( V * 0 > + B m 4 V ]di e 


[I ^ J ■ L * fc*. 


*IC ( 'l 


x o ,ldi « 



‘ J [ it x ic^ s' ‘ i S' ( 'i * 


^ • L * 


T 

y r y \ t. 

) x~, X 1 


T 

*1 _1 _ T 

X~,X 1 


dx 
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APHiNDIX G 


LINEARIZATION OF FORWARD FLIGHT EQUATIONS 
The nonlinear flap-lag-torsion equations of motion, Equation 

(U.22 ), 


G - [n(q)]{ \ } * [d( a ,q)][q} + [k(q)](q] + [t] - 0 (G.l) 

are linearized using a Taylor series expansion. The nonlinear mass, 
damping, and stiffness matrices and the forcing vector are defined 
below in indlclal notation. The range of the indices is 1 through 
M, where M is the total number of modes used in the analysis. 

[m(q)] - m^ + (G.2 ) 

[d(q,4)] -c lij+ dg^ ♦ d nj 

+ ^ C 2ijk + ^ijk" c Axijk + d 2ijk + d 3ijk/ q ^ q k 
+ (c 21ijk + ^ijk * d 32 ijk/‘ q / )q k (G,:5) 

(k(q)l - b UJ ♦ t lu ♦ k^ + a nj 

* ^ b 2ijk + k 2ijk + a 2ijk + a 3ijk/ ’ q /\ ^ G *^ 

[f] » f n + f A1 (G.5) 
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where 



- 

(A) 1 ty* 

[A) 

“2ijk 


- (A) 1 tlj I s 

(A) 

*2ijk 

\ ' 

■ :ai t d^j s 

t C% 1 ) S )(A) 

C A..iJk 


■ fA) T ([C A ^) S 

♦ Ccj jt ) S )[A) 

b lU 

X 

■ tA] T ((B 1 l S . 

■ [B C ] S )[A] 

f Ii 

n 

' (A) 1 (Fj) S 


f Ai 

- 

[a) t (f a } s 

• 


Ail other quantities in Equatione (0.2) - (0.5) are obtained directly 
fron the corresponding (upper case) system matrices, l.e., 

C 11J ' IAiI tc i) S (A) , 


etc. 

The derivatives of G, Equation (k.Zk) and Equation (U.3J*), 
are then defined as follows : 
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^ ’ b lin + t ll n + hin + a ii n 

+ (b 21nk + ^2 Ink * ^ink + a 3lnk/ q /\ 

+ (b 2ijn ^ijn* « 2 1Jn > 

* a 3ijkn q k q j 

(c 2ijn * ^ijn • c Ajcij n + ^ijn + d 3lJn/ q, 

+ a 3ijkn q k + d 32ijkn V q j 
+ “^ijn^ * (g.6) 

\ m 

^ C lin + d Sin + d iin 

+ (c 2lnk + *2 ink * C Axink + ^ink + d 3ink/ q /\ 
(c 21ink + ^2 ink + d 32inki q p q k 
(c 21ij n + ^ijn + d 32ijni q p q j » (G.7) 

5 . 

^ “lin ^ink q k * (g.6 ) 
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